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Abstract 


This thesis introduces a new specification and verification approach for dynamic systems. 
The introduced approach is able to provide type II error free results by definition, i.e. there 
are no hidden faults in the verification result. The approach is thus suitable to provide a 
reliable verification of safety critical systems. 


A new notion of set based consistency for dynamic systems with a given specification is 
presented. Therefore Kaucher interval arithmetic is used to enclose the measurement data 
in a bounded error sense. The resulting method is able to verify the specified behavior of 
a dynamic system against its measurement data even in the presence of noise and sensor 
uncertainty. Consistency is defined using the Kaucher arithmetic united solution set which 
leads to mathematically guaranteed results. 


It is proven mathematically that the desired property holds for a wide class of systems, in- 
cluding time invariant, interval type and hybrid systems, which can be used to describe even 
nonlinearities. Several extensions are introduced, leading to a new iterative identification and 
segmentation algorithm for hybrid systems which is able to handle even unknown switching 
times. In case the calculations can be done fast enough, the developed approach can also be 
used for the diagnosis of dynamic systems. 


The presented methods are successfully applied to several example systems, including theo- 
retic settings and a variation of different tank settings. 


The new theories, methods and algorithms developed in this thesis form the foundation for 
reliable safety analysis of highly automated safety critical systems. 


Zusammenfassung 


Diese Arbeit beschreibt einen neuen Spezifikations- und Verifikationsansatz für dynamische 
Systeme. Der neue Ansatz ermöglicht dabei Ergebnisse, die per Definition frei von Fehlern 
2. Art sind. Dies bedeutet, dass das Ergebnis der Verifikation keine versteckten Fehler enthal- 
ten kann. Somit können zuverlässige Ergebnisse für die Analyse von sicherheitskritischen 
Systemen generiert werden. 


Dazu wird ein neues Verständnis von mengenbasierter Konsistenz dynamischer Systeme mit 
einer gegebenen Spezifikation eingeführt. Dieses basiert auf der Verwendung von Kaucher 
Intervall Arithmetik zur Einschließung von Messdaten. Konsistenz wird anhand der vere- 
inigten Lösungsmenge der Kaucher Arithmetik definiert. Dies führt zu mathematisch 
garantierten Ergebnissen. Die resultierende Methode kann das spezifizierte Verhalten eines 
dynamischen Systems auch im Falle von Rauschen und Sensorungenauigkeiten anhand von 
Messdaten verifizieren. 


Die mathematische Beweisbarkeit der Konsistenz wird für eine große Klasse von Systemen 
gezeigt. Diese beinhalten zeitinvariante, intervallartige und hybride Systeme, wobei letztere 
auch zur Beschreibung von Nichtlinearitäten verwendet werden können. Darüber hinaus 
werden zahlreiche Erweiterungen dargestellt. Diese führen bis hin zu einem neuartigen iter- 
ativen Identifikations- und Segmentierungsverfahren für hybride Systeme. Dieses ermöglicht 
die Verfikation hybrider Systeme auch ohne Wissen über Schaltzeitpunkte. Die entwickel- 
ten Verfahren können darüber hinaus zur Diagnose von dynamischen Systemen verwendet 
werden, falls eine ausreichend schnelle Berechnung der Ergebnisse möglich ist. 


Die Verfahren werden erfolgreich auf eine beispielhafte Variation verschiedener Tanksysteme 
angewendet. 


Die neuen Theorien, Methoden und Algorithmen dieser Arbeit bilden die Grundlage für eine 
zuverlässige Analyse von hochautomatisierten sicherheitskritischen Systemen. 
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1 Introduction 


The fast technological development in computer engineering in recent years led to very pow- 
erful computing capacities that are now available at very low costs [Wil17]. As a result those 
chips are used in an increasing number of products to make them “smart” and to enhance user 
experience and functionality. These smart devices are interleaving the daily life of millions 
of people and are used for an increasing number of tasks [Gho17]. State of the art techniques 
are powerful and mature enough to take over even very complex and sensitive tasks - for 
example in autonomous cars, in flight assistance systems or in the control of critical infras- 
tructure. Tasks that can potentially harm human beings or destroy valuable infrastructure 
are called “safety critical” and special measures need to be taken during the development 
cycle to ensure correct operation of such safety critical systems [IEC10][ISO11]. 


These special measures are given by safety analysis methods. A very relevant property of all 
safety analysis methods is given by the amount of their type I and type II errors. Thereby 
type I errors (false alarms) denote the situation in which a safety analysis method evaluates a 
correct system to be faulty. The complementary condition is given by type II errors (hidden 
faults). In this case a safety analysis method evaluates a faulty system to be correct. Type II 
errors are of major importance in the context of safety critical systems. A faulty system that 
is evaluated to work correctly poses uncontrollable risk to the user and the environment. 
Thus there is a need for safety analysis methods that do not suffer from type II errors. 

In the context of this thesis, verification of dynamic systems means applying safety analysis 
methods in an offline setting to ensure consistency of the verification object (VO) with the 
specification. Guaranteed verification means that type II errors are impossible by design. In 
case the safety analysis is fast enough, it can be applied in an online setting which is then 
called “diagnosis”. Diagnosis can also be used to detect runtime errors. 

It is common opinion that there is currently no sufficient type II error free method available 
in the state of the art and the state of science [Kap16]. 


Currently safety analysis methods use mostly falsification approaches, e.g. methods from the 
field of testing. To achieve confidence about the absence of failures based on testing methods 
it is necessary to use a sufficiently large amount of test cases. This leads to the fact that safety 
analysis is more expensive as the development itself [Fos15] and costs are expected to rise 
further with increasing complexity of the tasks assigned to the technical system. 


2 1 Introduction 


Besides costs, current safety measures are often based on the experience of the responsible 
engineer or brute-force simulation approaches are applied [ZN09]. It is widely recognized 
that those methods will not keep up with the complexity given by currently developed or 
future systems [Ram17][Ott18]. 

The example which is predominant in public perception is located within the automotive in- 
dustry. Current systems like automatic cruise control or autonomous parking pilots are ana- 
lyzed by applying the previously mentioned experience and simulation based safety methods 
[Zan12]. Nevertheless it is known that those are not suitable for the arising challenges, e.g. 
in the context of autonomous driving [Wac17][Koe18]. 


A possibility to avoid this dilemma is to use formal methods that can proof specific properties 
of a verification object. A promising approach that gained great attention in recent years 
is given by interval arithmetic safety analysis methods, see among others [Uga03][Wol10] 
[San17]. Results obtained using those method are guaranteed to include all possible nominal 
system behavior as well as additional non-nominal behavior (so called spurious solutions). 
Due to this overapproximating property, there are no type I errors. However, for the same 
reason type II errors are possible by design. 


The goal of this thesis is to close the gap by developing a formal method for the safety analysis 
of dynamic systems that is guaranteed to be free of type II errors. 


2 State of Science 


There are numerous methods and approaches concerned with (safety) analysis of develop- 
ment results in different communities. Also, there is a broad terminology with respect to the 
verification question. The primary goal of this chapter is to build a basic conceptualization 
and the resulting terminology used in this thesis. This is necessary to follow the ideas and 
approaches introduced in later chapters. 

Furthermore, the most important and wide spread notions and methods used in engineering 
and engineering science are introduced and discussed. 


2.1 Conceptualization and Terminology 


System behavior analysis can be conducted with respect to different perspectives. This thesis 
addresses the verification of dynamic systems. A classification of other perspectives is given 
in Appendix A. To conduct the verification of dynamic systems it is necessary to define three 
components. First a description of the intended system behavior is set up. The next step is 
to define the concept of deviating behavior. Finally the developed real system behavior has 
to be assessed with respect to the intended behavior. 


2.1.1 Behavior Description 


The desire of the costumer needs to be documented in some kind of specification to allow 
any analysis in terms of verification. There are as many specification methods to define 
the nominal behavior as there are methods to check their fulfillment. The variety includes 
very formless approaches in human language [Mac95] as well as very formal definitions us- 
ing (runnable) models [A115] or special specification languages ([Par72][Spi89][Abr96]). The 
choice of a suitable specification formalism to be used in a project is a trade-off. The less 
formal a specification, the less effort is necessary to set it up, leaving the effort to the de- 
veloper who needs to interpret the specification. During the verification procedure it is nec- 
essary to interpret the specification which leads to a need for experienced experts [ZN09, 
p. 120][Raj13][Bal16]. The more formal a specification, the more effort is necessary to set 
it up. An advantage of formal specifications is that they force the specification engineer to 
capture the requirements in a precise and structured way. Standardized specification routines 
help to avoid careless mistakes during the setup [Par86][Hal90][Sch15a]. 

On the other hand all properties that should be covered need to be representable in the spec- 
ification, which can lead to requirements being impossible to be captured in a specific for- 
malism. However, due to the precision of very formal specifications it is possible to analyze 
them in a rather automated or “proof-like” way. 
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Throughout this thesis it is assumed that the specification itself is known, correct and rep- 
resents the whole functionality, behavior and all properties that are necessary to fulfill the 
customers desire. 


This work assumes a set based specification. It is assumed that it is possible to represent 
the desired behavior of a dynamic system in terms of a specific abstract set. This formal 
specification (FS) can be interpreted to include all dynamic system parametrizations that are 
able to create the intended behavior. The set of intended or desired behavior can be given by 
differently shaped sets e.g. a circle or a square. In case there is no variation in the desired 
behavior, the set consists only of one parametrization which is given by a distinguished point. 
Different possible specification sets are depicted in Fig. 2.1. 


x © 7 


Figure 2.1: Exemplary set based specifications (Point, Circle, Square) 


2.1.2 Behavior Deviation 


Implemented systems can show behavior deviating from the desired behavior due to several 
causes. From a very basic point of view it is possible to differentiate between mistake by 
misfortune, mistake by accident without intention and mistake by deliberate wrong-doing 
by an individual! 

The setting of this work tackles the second kind, mistake by accident without intention that 
can happen at every point during the development process. A wide range of expressions is 
used to differentiate the field of unintended behavior or unintended properties. However, 
different fields of research and profession are using different naming conventions. 

The naming convention used in this thesis is given in Fig. 2.2. 


Misfortune Accident | Wrong-doing 
n ds 


Figure 2.2: Failure terminology 


! This categories are inspired by Aristotle (384 - 322 BC) who thought about ethics and mistakes of human behavior 


[Res07]. 
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The foundation of all unintended behavior is given by the basic mistakes. The next instance is 
called fault and denotes the deviation of at least one system value from its intended value.This 
deviation can happen due to all three of the basic mistakes. If a fault leads to unintended 
system behavior it is called error. 

A system perturbed by a fault and a resulting error can still be able to operate correctly. Only 
if there is a persistent interruption of correct behavior the system is called to show a failure. 
This failure introduces a hazard into the environment the system operates in. The hazard 
can lead to consequences in the environment that potentially harm objects or even human 
beeings. 


Complementary, there is the concept of disturbance in a control engineering sense. The pro- 
cess of capturing real world data and transferring them into any control system is always 
superimposed by a process that creates a deviation between the real values and the measure- 
ment values [Fral6, p. 65]. This deviation is called disturbance or noise and every system 
needs to be adapted to the specific noise present in itself as well as in the particular environ- 
ment. 


2.1.3 Behavior Assessment 


The assessment of the verification object is done with respect to its behavior. Therefore it 
is necessary to set up the formal specification (FS) and additionally describe the behavior of 
the verification object (VO) using the same formalism. Both descriptions are assumed to be 
represented by a convex set. The notion of set based basic consistency that is used throughout 
this thesis is given in Definition 2.1. 


Definition 2.1 (Set Based Basic Consistency) 


A set based verification object VO is called basic consistent with its set based formal spec- 
ification FS if and only if there is an intersection between the formal specification and the 
verification object behavior: 


(FSnVO £z 0) & Basic Consistency. 


A special case is given by full consistency, which means that all behavior given by the formal 
specification is available in the verification object. 


Definition 2.2 (Set Based Full Consistency) 


A set based verification object VO is called full consistent with its set based formal specifica- 
tion FS if and only if the formal specification behavior is an subset of the verification object 
behavior: 


(FS C VO) & Full Consistency. 
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The inverse is given by inconsistency according to Definition 2.3. 


Definition 2.3 (Set Based Inconsistency) 


A set based verification object VO is called inconsistent with its set based formal specification 
FS if and only if there is no intersection between the verification object behavior and the 
formal specification behavior: 


(FSnVO = 0) & Inconsistency. 


This means that none of the available VO behavior is given in the specification. 


The resulting situations are depicted in Fig. 2.3. It is assumed that the formal specification 
FS is given by the blue circle. The set of real VO behavior is given by the green and red 
circles. Definition 2.1 and Definition 2.2 are fulfilled in the left and middle subfigures, leading 
to the verdict Basic Consistency and Full Consistency for the depicted VO and the FS. 
Definition 2.3 is fulfilled in the right subfigure, leading to the verdict Inconsistency for the 
depicted VO and the FS. 


Formal Specification FS Formal Specification FS Formal Specification FS 


o 


VO Behavior VO Behavior 


(Full Consistency) 


VO Behavior 


(Basic Consistency) (Inconsistency) 


Figure 2.3: Basic consistent, full consistent and inconsistent result of set based verification 


In the context of this thesis the verification object is considered to be correct if there is spec- 
ified behavior within the VO behavior. This is called "consistent behavior". 

For the ease of notation, the term Consistency is used as soon as there is “consistent behav- 
ior”, either due to Basic Consistency or due to the even stricter Full Consistency. The 
following considerations apply equivalently to both definitions. 


In general, the VO behavior is not directly available and thus needs to be captured by an 
approximation. If the behavior is given in terms of dynamic system parameters, the approx- 
imation can be calculated using identification methods [Lju99]. Therefore it is necessary to 
interact with the real VO to determine the underlying behavior. Assumption 2.1 has to hold 
to allow a successful identification. 


Assumption 2.1 (Persistent Excitation of the VO) 


The verification object VO is sufficiently excited to show all relevant behavior. 
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Only behavior of the VO that is triggered or excited is included in the approximation and 
can thus be analyzed [Ast95, p. 41][Ise10, p. 250]. Throughout this thesis it is assumed that 
Assumption 2.1 holds. 


There are two main set based calculation paradigms that can be used to determine the ap- 
proximation of the VO: underapproximation (^) and overapproximation (*). 

In case of overapproximation, there is spurious behavior in the resulting outer enclosure (see 
rectangle in the left of Fig. 2.4). If underapproximation is used, some VO behavior is missing 
in the inner enclosure (see rectangle in the right of Fig. 2.4). 
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Overapproximated 
Behavior VO* 


Figure 2.4: Set based overapproximation and set based underapproximation 


As the true VO behavior is not available, the approximated behavior is used to reason about 
Consistency. In case an overapproximation of the verification object behavior (VOT) is used, 
this leads to 


(FSNVO*F z 0) & Consistency”. (2.1) 


It is very important to note that (2.1) yields the verdict Consistency* that holds only for the 
overapproximated behavior VO *. It does not have to be valid for the true VO behavior. To 
relate the verdict with the true VO behavior, the overapproximating property 


(VO c Vo*) — Consistency C Consistency* (2.2) 
has to be taken into account, leading to 


(FS nVo*z 0) Consistency? = Consistency & (FS nVoz Ø). (2.3) 


It can be seen from (2.3) that the Consistency” verdict can not be extended to the true system 
behavior VO. This situation is depicted in the lower left field of Tab. 2.1. This property holds 
for both, Full Consistency and Basic Consistency. 


8 2 State of Science 


If an overapproximating method yields the result Consistency* and this result is general- 
ized to Consistency, it is possible that there is a “hidden fault” present in the system. This 
situation is called type II error and it is a severe problem in the field of safety critical systems. 
Type II errors are likely to harm people or the environment as the VO is showing wrong 
behavior but the supervising system assumes correct functionality. Corrective actions that 
are designed to prevent fault induced damage are not activated in case of a hidden fault. This 
type of fault can thus proceed and potentially harm people. In safety critical systems, type II 
errors need to be avoided under all circumstances. 


Therefore it is beneficial to use the underapproximation of the verification object behavior 
(VO), leading to 


(VO- C VO) => Consistency” C Consistency. (2.4) 
The resulting verdict Consistency” can be extended to the true system property: 


(FS QVO Æ 0) = Consistency” => Consistency = (FS NVO Æ 0). (2.5) 


If the underapproximation of the verification object VO” is consistent with the specification, 
it is guaranteed that the true verification object is also consistent with the specification (see 
Tab. 2.1, lower right field). Thus the verdict is guaranteed to be free of type II errors. Again 
this property holds for Full Consistency and Basic Consistency. 


The new verification approach developed in this thesis is based on the concept of underap- 
proximation to avoid type II errors by design. This essential property is necessary to solve 
the currently unsolved reliable verification problem of safety critical systems. 

However, the property comes at the costs of possible false alarms (see Tab. 2.1, upper right 
field). Instead of additional spurious solutions there are missing solutions generated by the 
underapproximation VO”. Even though there is consistent behavior in the VO, this behavior 
is not included in the underapproximation, leading to a false alarm (type I error). 

The situations in the remaining upper left field of Tab. 2.1 depicts the correct verdict 
Inconsistency than can be obtained using both, under- or overapproximation. This is due 
to the fact that there is no consistent behavior for the real VO as well as for both approxima- 
tions. 
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Verification Result: 
“Inconsistency” 


Table 2.1: Error types 


System is Faulty 


Nominal Behavior 
(Inconsistent) 


Correct Result 


System is Correct 


Nominal Behavior 
False Alarm 


nderapprox. 
O Behavior 
Loe 


Type I Error 


Verification Result: 
“Consistency” 


Nominal Behavior 
Hidden Fault 


A 
7 
2 


Z 


> | 


VO Behavior 


Type II Error 


Nominal Behavior 
(Consistent) 


O 
VO Behavior 


nderapprox. 
O Behavior 


Correct Result 


10 2 State of Science 


2.2 Interval Arithmetic Methods 


The overapproximating property introduced in the previous section can be achieved using 
the notion of interval arithmetic. Interval arithmetic is thereby used to enclose the effects of 
noise and epistemic lack of knowledge that is always present in real systems. A fixed lower 
and upper bound is used to describe a set of possible true values £irue that are associated 
with a given measurement value £meas: 


Ttrue = [Emeus m ô, Tmeas + ô]. (2.6) 


The maximum tolerance ô is the only parameter needed to set up the interval. This interval 
arithmetic notation of maximum deviation is widely used by sensor manufacturers [Kre95] 
and in the fault detection community (e.g. in [Arm09][Zai14][AI17]). The true measurement 
value is guaranteed to be included in the interval and it is guaranteed that the true value is 
never outside the interval. 

When interval enclosure is used on the measurement data, all succeeding calculations have 
to apply the notions of interval arithmetic to preserve the guaranteed properties. The basic 
property of interval arithmetic calculations is that all possible solutions are included in the 
result. Therefore, interval arithmetic results are able to create the introduced overapproxi- 
mation properties. The interval arithmetic solution set consists of real solutions and spurious 
solutions. Spurious solutions denote solutions that do not exist in the real system but are in- 
evitable artifacts that are introduced by interval arithmetic calculations and the final interval 
arithmetic (and thus axis parallel) enclosure of the real solution [Bau87]. 


Interval arithmetic is widely used for verification [Bal16] and diagnosis methods as type I 
errors (false alarm) are prevented by definition. Therefore, models of the nominal system are 
used to calculate a set of predicted outputs for the measured inputs of the system [Ven15]. 
This can be done by using intervals on the system parameters to calculate an interval range 
of outputs [Pui06][Mes10][Wol10]. The system is assumed to be correct as long as the mea- 
sured output values are within the predicted output interval, i.e. within the so called direct 
image. Due to the used outer enclosure of the prediction, type II errors are possible using 
this class of methods. 

An alternative approach uses the so called inverse image [Pui06] or feasible set [Cas14], also 
known as set-membership approach [Ing09]. In this case the input-output measurement data 
is used to calculate the set of parameters that is able to generate the observed mapping. This 
is possible if the system is linear or nonlinear but linear with respect to the parameters. If 
there is a member of the set of nominal models within the feasible set of the measurement, 
the measured data can be explained by the nominal model. This class of methods utilizes in- 
terval arithmetic identification based on outer enclosures. Therefore there are type II errors 
possible by definition. 

The feasible set resembles the solution set of the identification problem given by the mea- 
surement data that can be computationally hard to calculate [Hor13]. 

A wide spread possibility to approximate the solution set is given by subpavings using the 
SIVIA algorithm (see [Jau01, p. 45ff][Pui06][Mes10]). The bisection approach of SIVIA leads 
to a large set of different intervals with various size. Even though this result is very precise, 
the great amount of intervals leads to a complex handling and to long calculation times. 
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A more efficient approach is given using a zonotope? representation as shown in [Ing09]. 
The question how this approximation can be calculated in an efficient way is still an active 
research topic. Latest results [Koc19] use sparse polynomial zonotopes. Additionally, there 
are methods to reduce the solution set by pruning spurious parts if possible. For example 
the approach presented in [Wol10] uses measurement data to reduce the overapproximated 
state set. Nevertheless, all approximation methods use outer enclosures and are thus prone 
to type II errors. 

Besides the field of diagnosis, different interval arithmetic approaches are used for state es- 
timation ([Jau01][Ram09][Mes10][Efi13] [Kre16][Kre18] [Wan18]) and control [Rau06]. None 
of these methods addresses type II errors. 

Therefore it is necessary to develop a new verification method that is able to guarantee the 
absence of type II errors. This can be achieved by calculating an inner enclosure that consists 
of a subset of the real solution as shown in the previous section. This thesis utilizes a special 
extension of interval arithmetic called Kaucher interval arithmetic that provides powerful 
theories to calculate the necessary inner enclosures. 


2.3 Governing Complexity: Time Variant and Hybrid 
Verification Approaches 


Linear time invariant (LTI) systems form the base for the most system theoretic methods and 
approaches. Nevertheless, real world systems are normally neither linear nor time invariant. 
It is thus the question how to handle the complexity inherently given in real world systems. 
Some nonlinear systems can be handled using nonlinear theory which provides methods 
for analysis and control [Kha15]. However those approaches are often subject to very strict 
preliminaries and only applicable for a narrow class of systems. Therefore linearization is 
often used in practice to handle nonlinearities. Another possibility is to model nonlinearities 
by using time variant parameters in a linear system [Ble11]. In this case, the set of feasible 
parameters can be bounded using interval arithmetic or any other set definition. It is also 
possible to split the nonlinear dynamic into a sequence of linear dynamics [Oza14] that are 
activated by a superimposed switching mechanism. The resulting model belongs to the class 
of hybrid systems. Hybrid systems can also be used to model time variant linear systems 
with piecewise constant parameters. Therefore the different piecewise constant parameters 
are represented using an individual dynamic subsystem each. 


There is a wide theoretical framework available for hybrid systems [Eng02][Mah10]. The 
topic of hybrid verification is subject to current research in the cyber physical systems com- 
munity (see e.g. [Sch15b][Kap16][Sch17a][Ara17] [Bar18][Har18][Lau18]). There are also 
large research clusters in this area [AVA19][ENA19]. 

To apply the set based approach introduced in this thesis in a hybrid setting, it is necessary to 
use hybrid identification methods to determine the unknown system parameters from given 
measurement data. The most relevant hybrid identification approaches given in the literature 
are introduced in the following. 


? A zonotope is a convex polytope that is point symmetric with respect to its center. 
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The algebraic approach provided in [Vid08] interprets the identification setting as a geomet- 
ric problem. The measurement data as well as the parameters are interpreted as vectors that 
have to be perpendicular in case the identified parameters match the true parameters. The 
goal is to find the parameter vector with minimal projection on the family of all measurement 
vectors. A Bayesian approach based on stochastic properties is given by [Jul05]. The number 
of models and the model order need to be known a priori for this procedure. A cluster based 
approach was developed by [FT03] were machine learning methods are used to form groups 
of similar behavior. Identification methods based on optimization were developed and pre- 
sented in [Mün05][Bor09][Lau18]. The bounded error approach introduced in [Bem05] and 
used for time variant systems in [Bra16] assumes errors that are characterized by their max- 
imum value. Even though this is close to the basic interpretation used in this thesis, those 
approaches do not use interval arithmetic notations. Therefore they lack the guarantees that 
are necessary in the safety critical context of this thesis. 

A greedy approach based on [Oza12] was developed in [Die13a] and [Die13b]. This approach 
is different from the others as it is the only one that uses a multi-step prediction error in- 
stead of the common one step prediction error. It was extended to cyber-physical systems in 
[Sch17a] and to Kaucher interval arithmetic in [Sch19]. 


2.4 Other Common Verification and Falsification 
Approaches 


There is a wide range of verification methods with different degree of abstraction and for- 
malization. 

A strong diagnosis community is active in the control engineering field (see e.g. [Ise93] 
[Sch03][Ven03a][Ven03b][Ven03c][Bla06][Arm09][Sch09][Pui10][Ise11][Cac13][Zol14]). 
Also, a large testing and verification community formed in the information technology com- 
munity. This leads to a great range of verification and validation methods, unfortunately 
using a similar terminology (e.g. [Bar78][Boe84][Hay86] [TF91][Bar05]). Three of the most 
relevant approaches are introduced in this section. 


2.4.1 Testing 


A wide spread - if not the most wide spread - approach to assess the properties of a technical 
system is given by testing. Testing is a classic falsification method, aiming on the detection of 
counter examples that do not show the intended behavior. It is therefore necessary to define 
a well-chosen set of test cases, consisting of the system state at the beginning of the test case, 
inputs that are applied during the run ofthe test case and outputs that are expected to appear 
during or at the end of the test case. If the system under test (SUT) shows outputs deviating 
from the expected outputs, the test case is called a "failed test" and further inspections of the 
test case are necessary to identify the reason. In contrast to verification methods, each result 
is only valid for the specific applied test case. It is long known that the confidence of the 
results can only be increased by increasing the amount of test cases [Fut89, p. 3]. 
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2.4 Other Common Verification and Falsification Approaches 


Testing can be applied at different stages of the development cycle and at different abstraction 
levels. The state of the art testing scheme is given by the V-model (see e.g. [Web09][ZN09] 


[Raj13][Ott18]) as depicted in Fig. 2.5. 
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Figure 2.5: V-model diagram (based on [Ott18]) 


The left part of the V-model is the specification branch. It is applied in a top-down fashion. 
The specification is initialized on the highest level representing the customers desire. Then it 
is propagated to lower levels and refined to match the degree of formalization of each level. 
This structured procedure includes the decomposition of the overall task in several sub tasks 
accompanied by the definition of interfaces between the sub tasks. 
After every refinement step the result is checked against the superimposed level to verify 
that both requirements are consistent. The final (software) specification is implemented at 
the bottom level. The resulting SUT is checked against its specification at the same level. If 
the test succeeds, the function is used in higher levels and combined with other parts to meet 
superimposed specifications. If a test fails, the system under test is rejected to the previous 
level. The complexity of the SUT increases with rising level on the right verification branch 
of the V-model. It is also possible that a failed test at high levels (e.g. at system or acceptance 
level) leads to changes in the respective high level specification. This results in the need of 
repeating the whole development process starting with the changed high level specification 


[Raj13][ZN09]. 
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There are several technologies that are used in different phases of the V-model, sometimes 
leading to additional branches. Fig. 2.5 includes cutting edge technologies like rapid proto- 
typing, component HiL and cluster HiL that are used to speed up the time effort to perform 
a complete cycle of the V-model. 

An important part of the testing process is the test case generation. There are several possi- 
bilities to set up the test cases. If there is experience with this kind of product, it is likely that 
an existing test case database can be reused and adapted [Sax08][Bal16]. Another straight 
forward approach is to determine the feasible range of all input variables and divide this 
range in so called “equivalence classes” [Utt06]. The test cases are then formed by combining 
one representative (or the minimum and maximum values) from each equivalence class with 
all possible representatives of the other input variables [Utt06]. If it is not possible to form 
equivalence classes, it is also possible to sample the valid range of a variable randomly or in 
equidistant steps. This approach is straight forward and easy to understand but it will lead 
to very large sets of test cases with increasing number of variables in the system or if a fine 
resolution of the variable ranges is necessary (i.e. [Hei05][AI15]). 

The runtime for testing or simulation rises with the number of test cases, leading to large 
computation times. If equidistant sampling is used, it is possible that much calculation time 
is spent assessing “uninteresting” regions of the input space or that “interesting” regions are 
only covered with few test cases. One major drawback of the testing approach is that all test 
cases need to be redone if the SUT changes. As there are frequent iterative changes during 
a development cycle (i.e. more than one iteration of the V-model is necessary), the test cases 
need to be applied several times which leads to even longer computation times. 


2.4.2 Reachability Analysis 


By including more system theoretic knowledge, testing can be developed to reachability anal- 
ysis. The basic requirement for this purpose is a specification that includes some kind of state 
space the system is operating in. Further, the specification has to define forbidden areas in 
this state space. 

The main idea of reachability analysis is to determine whether there is a sequence of inputs 
that leads the SUT to enter the forbidden area. If an input sequence leading to a forbidden 
area is found, the SUT is falsified and needs to be improved to match the requirements (see 
among others [Bha04][Alu06][Mit07][Alt08][Don10]). Reachability analysis is also aiming 
on finding a counter example which has the same basic problem as the testing scenario: no 
detected counter example does not mean there is no fault present in the system as faults can 
be hidden in uncovered parts of the state space. 

Runtime limitations restrict all methods to a finite number of samplings and thus to partial 
coverage ofthe state space. There are smart coverage criteria available that allow an effective 
calculation of the most important regions of the state space e.g. using rapidly exploring ran- 
dom trees (RRT) (see i.e. [Bha04][Kap16][Pan17]). These approaches can be extended to the 
so called method of star discrepancy [Dan11] or by applying the underminer method [Bal16]. 
Nevertheless, reachability analysis is still a falsification method, based on the specification of 
the faulty case (forbidden areas). Therefore reachability analysis is not suitable to solve the 
verification problem addressed in this thesis. 
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2.4.3 Formal Verification 


One possibility to avoid the counter example problem is given by more abstract approaches 
from the verification field. Those formal methods are used to reason about system proper- 
ties in a mathematical rigorous way. To apply formal methods, the VO needs to be trans- 
ferred to a strict mathematical notation, e.g. by using the Z specification language (see e.g. 
[Spi89][Bro05, p. 325]) or Prolog [Bro05, p. 334]. The formalized VO can then be used to carry 
out mathematical proofs showing that specific system properties hold in all operating condi- 
tions. The proof is thereby conducted by a so called theorem prover. 

This approach is very powerful as the results are valuable and mathematically sound. Never- 
theless, formal proofs can only be done for very distinct properties. Furthermore the methods 
need very long runtime even for “small” academic problems. This leads to still unsolved run- 
time issues for real world problems ([Bro05, p. 325][Bar18]). 

A basic problem that cannot be omitted is that formal proofs cannot be conducted on the VO 
directly. Therefore the results hold only for the image of the VO which is given in the used 
formalism. Mistakes that are introduced when the system is transferred from the real world 
into the modeling formalism cannot be detected. 


2.5 Scientific Gap and Related Research Question 


Testing is the state of the art for current systems and is successfully applied in various com- 
munities. Nevertheless there are current systems, e.g. autonomous driving functions, that 
show a number of relevant scenarios that cannot be covered by testing or simulations. Even 
if this was possible, falsification methods cannot prove the absence of all faults. They need 
to stop at some point and have to assume that no undiscovered fault is present in the system. 
There are established methods available that use interval enclosures to mathematically bound 
the system behavior. Those methods use classic interval arithmetic, leading to overapproxi- 
mating properties. Overapproximating methods are able to provide type I error free results, 
meaning that there are no false alarms generated by the method. Nevertheless, the overap- 
proximation can cover missing behavior, leading to an undetected hidden fault. In the case of 
safety critical systems, hidden faults (type II errors) can lead to severe consequences threat- 
ening human life. It can be concluded that the verification of safety critical dynamic systems 
is currently not solved. 


A new verification method has to be developed to close this gap. This method has to be free 
of hidden failures, meaning that there are no type II errors. Therefore the specification is 
assumed to be formally given in terms of a set of dynamic system parameters. The behav- 
ior of the VO has to be given in the same formalism, leading to an identification problem. 
Safety critical systems are often implemented as embedded systems that consist of a closely 
connected discrete event system (the controller) and a dynamic system (the plant). Thus it is 
necessary to develop a hybrid identification method that is able to provide the desired guar- 
antees. 

The comprehensive research question tackled by this thesis is: 


“How can the consistency of highly automated safety critical dynamic systems be evaluated 
by a guaranteed verification method?” 


3 Methodical Approach and Mathematical 
Preliminaries 


Considering the state of science as well as the current and future challenges of system the- 
ory there is a need for a new verification method. The rising importance of safety critical 
systems emphasizes the need for formal methods that target type II errors. This thesis intro- 
duces such a formal method based on the notions of interval arithmetic, extended to Kaucher 
interval arithmetic. First the necessary notations and definitions are given to provide a sound 
theoretical base for further considerations. 


3.1 Mathematical Preliminaries 


All methods introduced in this thesis are based on interval arithmetic, appended by the prop- 
erties of Kaucher interval arithmetic. In the following, the basic properties and notations of 
interval arithmetic are introduced. 

The goal of this chapter is to provide a brief overview of interval arithmetic that is neces- 
sary for this thesis. The interested reader is referred to [Bau87][Rze08][Roh12][Sail4] for an 
extensive coverage of the topic. Throughout this thesis the well known notation of interval 
arithmetic extended by Kaucher interval arithmetic introduced in [Kup95][Sha96] is used. 


3.1.1 Basic Interval Arithmetic 


Interval arithmetic was initially developed to handle numerical calculation errors due to float- 
ing point calculation used in computer algebra systems [Apo67]. It gained popularity outside 
the numerical community with the rise of electronic computing in various fields. When mea- 
surement data is used in computing - as it is normally the case in engineering and natural 
science - faults are already created by the measurement process itself [Kre95]. Furthermore, 
the used values are given as samples at discrete time steps k. Every measurement Ymeas,k is 
compromised by some noise €x that leads to a deviation between the real value Ytrue,k and 
its measurement 


Ymeas,k = Ytrue,k + Ek. (3.1) 


It is possible to define intervals around the measurement that are guaranteed to include the 
real system value 


Utrue,k € [Ymeas,k = à, Ymeas,k = ô] (3.2) 


if the maximum 6 of the absolute noise is known i.e. Vk : |e;| < 6. 
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Suitable values of 6 can be determined from the data sheets provided by the sensor manufac- 
turers. 


The definition of a classical interval type variable x as given in [Sail4] is 
zx := |z, 7| = {rER|z <x <7}. (3.3) 


This definition includes all real numbers that are between or on the infimum x and the supre- 
mum z. In case x > 0 and 7 > 0 the interval is called positive interval. If the infimum is 
negative (x < 0) and the supremum is positiv (x > 0) i.e. if 0 € x, the interval is called zero 
interval. A negativ interval is given if x < 0 and 7 < 0. One last definition covers the case 
of supremum and infimum being the same, i.e. x = 7, which is called a degenerated interval 
[Dja17] or point real interval [Sai14]. 

The set of so called proper intervals is given by 


IR := {æ = |z, 7| |x < Zand z,7 € R}. (3.4) 


Despite this infimum-supremum notation, each proper interval can be given using the cen- 
ter 


(z+ x) (3.5) 
and the radius 
SA: L(rÉr-a). (3.6) 
of an interval. The interval can now also be stated in the center-radius notation 
== (xe mA). (3.7) 
Furthermore, it is important to introduce the interval type vector matrix notation based on 
[Jau01]. Vectors and matrices are written as capital letters X and interval values are given 


in bold font x, leading to interval matrices denoted as X. An interval vector X is defined as 
cartesian product of n closed intervals that includes a subset of the real numbers R: 


X — aU x x x... x 2), with) = lz, au fori € {1,2,...,n}. (3.8) 


This notation can be interpreted as projection of the i-th interval component a“) to the i-th 
axis of the vector space. An illustration for n = 2 and n = 3 is given in Fig. 3.1. 
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2) 
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ET 


Figure 3.1: Examples of the graphical representation of X € IR? (left) and X € IIR? (right) 


An (m x n), m,n € N, interval matrix A can be interpreted as subspace of R”*”. Again it 
is defined using the cartesian product of m - n closed intervals: 


go s 40 
A= : : (3.9) 
amd) ... afm) 
=a) x at) x... x gi) 
E (a) (3.10) 
with 1 <i € m,1 < j € n. The center matrix is defined element-wise as in [Hla14] to 
T dorus T 
A, ern (afe) =} (60 + a) ei 
as well as the radius matrix 
ij 1 Ps = 
Ay €R™™*” ; Ou) x (ates) — a) (3.12) 


with 1<i<m1l<j<n. 


The four basic arithmetic operations addition, subtraction, multiplication and division, i.e. 
* € {+,—,-,/}, are well defined for intervals. The general application of each operator on 
two interval values x = [r, 7] and y = [y, y] is given by 


Bea SHS ey | eee oy sys go} (3.13) 


according to [Apo67], which leads to the interval type calculation rule 


[x, 7] x [y, y] = [min (£ xy, z*y, Try, 2*7), 


(3.14) 
max (zxy, zy, Try, 2x9). 


The different elements within the min(-) and max(-) operations are due to the fact that the 
combination of all extreme values need to be taken into account. 
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Unfortunately, this property also causes two major drawbacks of interval arithmetic, the 
dependency effect and the wrapping effect. Those effects are explained in Example 3.1 and 
Example 3.2. 


With the assumption 0 ¢ [y, y] it is possible to explicitly state the four basic operations as 
in [Apo67]: 


lz, z] + [y, v] = [e +y, z +7] 
lz, z] - [y, y] = [z- 9, 2 - y] 
[x, x] - [y y] = [min (xy, xy, Ty, Ty), max (xy, 2J, TY, Ty) | (3.15) 


It is shown in [Apo67] that associative and commutative property hold for interval values as 
well. However, the distributive law is not applicable anymore and needs to be changed to 


a-(y+z)Cau-ytau-z (3.16) 


which is known as the subdistributive property for the interval values x, y and z. 


While evaluating an expression, every appearance of an interval variable is treated individ- 
ually as if it was independent from its other occurrences. Multiple occurrences of the same 
variable thus lead to a widening of the enclosure. This property is called dependency ef- 
fect and is illustrated in Example 3.1. One approach to mitigate the dependency effect is to 
reformulate the expression such that each variable occurs only once, if possible. 
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Example 3.1: 


Assume the function 
1 
I«)=1+ = (3.17) 


The resulting enclosure for the interval x = [1, 3] can be calculated straight forward using 
the interval arithmetic definitions of (3.15) to f(x) = 3; 2] , which matches the true range 
of the function within the interval. It is depicted by the blue dotted frame in Fig. 3.2. 

If (3.17) is reformulated such that there are multiple occurences of x e.g. 


~ r1 


(a) = 


(3.18) 
Hi 


the interval arithmetic evaluation yields f (a) = E 4]. It can be seen that this is a large 
overestimation of the true values of the function within x, depicted by the dashed frame in 
Fig. 3.2. 


Figure 3.2: Dependency effect based on f(x) and f(x) 
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Another effect occurring with interval calculations is the wrapping effect. This effect is 
caused by iterative calculations based on previous overestimations. Such iterative calcula- 
tions are e.g. necessary to solve an initial value problem or to evaluate a state space equation. 
An example for the wrapping effect is given by the initial value problem of Moore [Bau87] 
and is depicted in Example 3.2. 


Example 3.2: 


Assume the initial value problem for the differential equation 


Y (t) = e a) Y(t) (3.19) 


with Y (0) = [[—e, e], [1—e, 1+ ]]", e > 0. When the equation is evaluated, the result- 
ing solution set needs to be framed by axis parallel enclosures after each step. The solution 
sets fort; — i- At, At = %, i € {0,1,2,3} are depicted in Fig. 3.3. It can be seen, that 
the overestimation is continually increasing, as the inherited overestimation is passed on and 
used as base for further calculations. 


Figure 3.3: Example for the wrapping effect using the initial value problem of Moore [Bau87] 


The wrapping effect can reach a serious extent even after only one iteration. A minimal 
example illustrating the extent of the problem after two steps is given in Example 3.3. 
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Example 3.3: 


This example clarifies the effect of interval calculations in the setting of a proportional gain 
system with unknown gain p: 


up=y. (3.20) 


The system setup is depicted in Fig. 3.4. The input and output ranges of the system are given 
and can be included in the intervals u = [2, 3] and y = [4, 9]. The goal is to calculate the 
gain that maps all possible input values u € u to the specified output range y. 


p 
u= (2, 3] y = (4, 9] 


Figure 3.4: Proportional gain system with proper solution 


Using the introduced interval arithmetic calculations leads to 


p= Yu 


= (waz 2] 


_ [4 9 
[4 
e (1.3, 4.5]. (3.21) 


Re-substituting p into the system equation yields 


y-—pu 
= [1.3, 4.5] [2, 3] 
= [2.6, 13.5] 4 [4, 9] (3.22) 


which is a strong overestimation of the genuine output range. 


The example shows that the system parameter calculated from input and output ranges can- 
not be used to reason about the parameter set that is suitable to map the given input on the 
given output. The wrapping effect is caused by considering the combination of the extreme 
values of both intervals. Nevertheless, when regarding the task at hand in Example 3.3, the 
goal is not to find all possible gains connecting the two intervals but to find those gains 
reasonably connecting “the most” elements of the intervals. This slight but very important 
change is illustrated in Example 3.4. 
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Example 3.4: 


Assume the setting of Example 3.3. The parameter interval is calculated as before using 
p = y/u. The question is now how many pairs (u,y)|(u € u),(y € y) exist for each 
parameter p € p. Therefore the intervals u and y are divided into equidistant parts of Au = 
Ay = 0.0001. The resulting 10'001 discrete samples of us are combined with the resulting 
50'001 samples of ys to calculate the connecting parameter p, = ys/us. The histogram 
formed by 500'060'001 values of ps is depicted in Fig. 3.5. It can be seen that the extreme 
values of the outer enclosure of the solution p = [1.3, 4.5] are only connected by a single 
input-output pair each. On the other hand, there is a plateau between p; = |2, 3] that 
connects a nearly constant number of input-output pairs. Substituting this interval value 
into the system equation leads to 


yi = piu 
= [2, 3] [2, 3] 
= (4, 9 =y (3.23) 


which is exactly the given output range. The interval p; is an inner enclosure of the solution 
set of p — y/u. 


2.500,000 
2.000,000 
1,500,000 
1,000,000 
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1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 
Parameter p 


Figure 3.5: Distribution of parameters in the proper case 


The plateau in Fig 3.5 contains those parameters p that are able to map any u € u to a value 
y € y. Note that not necessarily all values y € y have to be met by pu. The contour of the 
histogram given in Fig. 3.5 can also be analytically calculated. The derivation of the exact 
distribution is given in Appendix B. 

The property leading to the wrapping effect displayed in Example 3.3 and Example 3.4 is 
the non-existence of an inverse element in classical interval arithmetic [Apo67].The inverse 
element in the real numbers is defined with respect to an operation and denotes an element 
that maps itself on the neutral element of this operation (see [Bro08, p. 340]). 


3.1 Mathematical Preliminaries 25 


Thereby the neutral element is also defined with respect to the same operation and denotes 

an element that maps each other element on itself (see [Bro08, p. 339]). 

There are neutral elements in classical interval arithmetic. For example the neutral element 

for addition is given by e, = [0, 0] and for multiplication by em = [1, 1]. Applying the 

neutral elements to an arbitrary interval value r = [r, 7] leads to 
r+e=[r+0,r+0]=[r, 7] (3.24) 
rem —[r- 1l, 7-1] =[r, T]. (3.25) 


However, in general there is no inverse element as can be seen in the following: 


r+(=r)= [r+ (7r), F+(-2)] # eaifr AT (3.26) 
p (+) = B | Bere (3.27) 


Equations (3.26) and (3.27) hold if and only if r = r which means that r is a degenerated 
point real interval [Apo67]. 


3.1.2 Kaucher Interval Arithmetic 


It is beneficial to define an extension to interval arithmetic that provides the existence of an 
inverse element for all arithmetic operations. By using Kaucher interval arithmetic [Kau80], 
the set of proper intervals can be extended by the introduction of a new set of so called 
improper intervals. These improper intervals are defined complementary to classical inter- 
vals: 


KR = {x = |x, z] |%<awandz,7e R}. (3.28) 


The set of all proper and improper intervals is given by IIR" = IIR U KR and is depicted in 
Fig. 3.6. The set of point real intervals is depicted as diagonal line in the figure. The set of 
proper intervals IR is formed by the half plain above the point real line. It can be seen, that 
the improper intervals KR complete the IR* by covering the half plain below the point real 
line. 


The definitions of the basic arithmetic operations x € {+,—,-, /} need to be adapted to hold 
as well for classical as for Kaucher interval arithmetic [Sha02]. 

Therefore, the definition of the negative part x” and the positive part x” of a real number x 
is given by 


xz? = max(—2,0) and xz? = max(z,0). (3.29) 


The four classic operations can thus be written as follows: 


x+y=]|z+y, 7+7] (3.30) 
z-y=|2-9T-y (3.31) 
x: y = [max (2®y®, 7°y°) — max ee) ; 

max (a Puy, gg?) — max D y* ja y y] (3.32) 


rz/y-—m: [1/9, 1/y| , for y: gy > 0. (3.33) 
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Figure 3.6: Geometric interpretation of IR*. The diagonal line represents point real values (based on [Sai14, p. 18]) 


In addition the two unary operators 


opp ([z, z]) = [-z, —z] (3.34) 
dual ([z, Z]) = [z, x]. (3.35) 


are introduced to toggle between proper and improper intervals. Using this operators leads 
to the definition of inverse elements in Kaucher interval arithmetic: 

a + opp (z) = |z, 2] + [-z, -=] = [0, 0] =: 0 (3.36) 

x/dual (x) = |z; z|] - [1/z, 1/z] = [1, 1] =: 1. (3.37) 


Those comply with the classic interval analysis definitions if all used intervals are proper 
[Sha02]. 


It is hard to imagine the nature of an improper interval as it is neither empty nor does it 
include the same values as a proper interval with inverse borders. A possibility to grasp an 
idea of the nature of an improper interval is given in Example 3.5. 
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Example 3.5: 


Assume the proportional gain setting of Fig. 3.7 which is similar to Example 3.3 but with a 
different output range y. 


p 
u = [2.3] y = |4, 5] 


Figure 3.7: Proportional gain system with improper solution 


The question is again which values can be used as gain p € p that maps all input values 
u = |2, 3] to the output range y = |4, 5]. The intervals are again divided into equidistant 
parts of Au = Ay = 0.0001. The resulting 10'001 discrete samples of us are combined with 
the 10'001 samples of y, to calculate the connecting parameter p, = ys/us. The resulting 
100'020'002 parameter values are used to set up the histogram given in Fig. 3.8. 
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Figure 3.8: Distribution of parameters in the improper case 


In this case it can be seen that there is no plateau in the histogram. However, there are two 
edges at pı = 5/3 and pa = 2. Substituting pı into the system equation leads to 


Yor = piu 
= 5/3 [2, 3] 
e [3.33, 5] Z [4, 5]. (3.38) 
The evaluation for p2 yields 
Up; = pau 
= 2[2, 3] 


~ [4, 6] 4 [4, 5]. (3.39) 
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It can be seen that both parameters are able to map some values of the input range u into 
the output range. Nevertheless there is not a single parameter that can map all values of u 
to y. This observation can be combined with the fact that the inner enclosure is an improper 
interval: 


p = y/dual(u) 


=wal E l 

[45 

[Ey 

x T9, 1.7] (3.40) 


Therefore improper intervals can be interpreted as solutions of a setting with “eroded plateau”. 


3.1.3 Interval Type Linear Equation Systems 


The introduced interval arithmetic considerations can now be extended to a vector ma- 
: i T T 

trix notation. Assume there are T measurement values (uj), , = [uı,ua,...,ur] and 

(Uk) La = [yı,ya,...,yr] , containing a valid range for each sample k € (1,2,..., T}. 

Each suitable parameter p € p has to comply with all input and all output ranges. This 

problem can be stated as an interval type linear equation system 


up = Yı 
up = % 
(3.41) 
upp = YT 
or more general, for vectorial input ES ee E , scalar output yi, and n parameters 
Tm 
[p 9), sez ,p(] : 
1 2 n n 
u pO + u D je dub + u pi ) — yi 
50 4 2) (2) n) (n) —_ 
up us p nus us p — y 
2 i Í ä (3.42) 


up) + Up p? LT ur p"? = yr. 
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The variables are used to set up the regressor matrix A € IIR T *"). the measurement vector 


B € IIR? and the respective parameter vector X € IR**! with 
A- (a?) = (u) (3.43) 

1«i«T, 1<j<n 1<k<T, 1<j<n 
Be pum = (Ye)ıcker “ey 
X= (29) = (n?) (3.45) 

1<j<n 1<j<n 
The interval type linear equation system can thus be stated as 

A-X=B. (3.46) 


This system can be interpreted as the collection of all point real linear equation systems that 
can be formed from the enclosed interval values [Sha96]. 

Equation systems with a regressor matrix of dimension T = n are called quadratic. Dimen- 
sion T « n stands for an underdetermined and dimension T' > n results in an overdeter- 
mined equation system. Underdetermined systems do not carry enough information to solve 
the problem unambigiously. This thesis focuses on overdetermined systems which is the 
most relevant case when regarding reasonable measurement times T and system orders n. 


The inverse of a quadratic point real matrix A is defined if the matrix is non-singular i.e. A^! 
exists if det(A) Z 0. Analogously a quadratic interval type matrix A is non-singular if all 
point real matrices contained in the interval matrix are non-singular ie. det(A) #0 VA € A 
[Sha14]. 

For overdetermined point real systems the criterion changes to a rank condition. The point 
real matrix A € R(T”) with T > n is said to have full rank if rank (A) = n. For interval 
type overdetermined systems, this condition again changes to rank (A) — n, VA € A. This 
means that all point real matrices included in the interval matrix need to show full rank. 
Determining the rank of an interval type matrix is in general an N P-Hard problem [Sha14]. 
Nevertheless, several criteria to check if an interval matrix has full rank were collected in [Sta16] 
based on [Sha14]. An introduction of the most relevant ones is given in Appendix C. 


The full rank condition is connected with persistent excitation according to Assumption 2.1. 
If the used input signal provides persistent excitation, the regressor matrix has full rank 
[Sha14][Lak14]. This leads to Assumption 3.1. 


Assumption 3.1 (Rank of the Regressor Matrix) 
The interval type regressor matrix A shows full rank according to [Sha96]. 


Throughout this thesis it is assumed that all specifications and measurements lead to an in- 
terval regressor matrix A that has full rank. 
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In general there is no unique, component wise point real solution vector X for such an in- 
terval linear equation system. Instead, (3.46) is solved by a set of point real solutions X`. The 
elements of the solution set X € Y? depend on the specific interpretation of (3.46). This 
interpretation is done by the interval quantors V and 3 as explained in [Sha02]. The notation 


of 


V[a, a] — 3 |b, b] (3.47) 


means that x has to solve (3.47) for all elements of (a € Rla € [a, a]) but only for at least 
one specific element of {b eRlbe [b, b] } [Sha02]. 

Each element of a vector or matrix can be assigned with an individual quantor, i.e. it is pos- 
sible to precisely define a specific solution set for the interval type matrix equations. Vectors 
and matrices with assigned quantors are denoted as B* or A*, respectively. A vector or ma- 
trix containing only the elements with assigned V quantor are denoted by BY and AY, the 
elements assigned with an 3 quantor are given by B? and A7. To split an assigned vector 
B* or a matrix A* depending on the quantors, the dualization of the intervals as given in 
(3.37) has to be used. 

According to [Sha02] the splitting is different for matrices and vectors and it is given by the 
following relation: 


Vector: B* :— dual (B") +B? (3.48) 
Matrix: A* := A” + dual (A7). (3.49) 
The specific elements of the matrices AY and A? are given by 
= €(1,3) fe = 
avis) = 1% mm (3.502) 
0 , else 
"m €(i,5) if ¢ =F 
axis) — 1^ iS (3.50b) 
0 , else. 
The elements for the vectors B" and B? are given by 
pe if € = 
po = | (3.51a) 
0 , else 
! bO ee 
PO — pe (3.51b) 
0 , else. 


The most general solution set is given by a mixed assignment of both quantors to the interval 
matrix A as well as to the interval vector B. The resulting AE-solution? set ^ 4 p according 
to [Hla14] is given by 

Dar (Af, BS) = (X €R™| (3.52) 
(VAY e AY) A (VBY e BY) A (34° € A?) A (AB? € B®) : ((AY + A?) X = BY + B?) } 


3 


The genuine notation of [Hla14, p. 2] is: 
X € R” is an AE-solution if VAY € AY, VBY € BY,34? € A?, 3B? € B? : (AV -A2)X = BY BP. 
This notation is slightly adapted for the sake of readability. 
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Based on the general AE-solution it is possible to define four distinct solution sets as given 
in [Sha96][Fie06][Hla14]. 


Definition 3.1 (United Solution Set 5 --) 


The united solution set is formed by all solutions of any of the point real systems A-X = B 
with A € A and B € B that are included in the interval system: 


S (A, B) := (X € R^| GA€ A) A (ABE B): (A. X =B)}. (3.53) 


Note that not all A € A can match any element B € B by multiplication with any X € » 744, 
and that not all B € B can be calculated using any X € » 744 and all available A € A. Fur- 
thermore, the set 5 = is not necessarily connected and not necessarily constrained by bor- 
ders parallel to the coordinate axes. Using an enclosing interval X 2 5 744 will likely create 
spurious solutions. 


Definition 3.2 (Tolerable Solution Set 5 ^...) 


The tolerable solution set includes all values of X that solve the interval type linear equation 
system regardless of the chosen point real matrix A € A. This means the solution holds for 
all included point real matrices: 


YO (A B) = {X eR"|(WAe A) ^ (IB € B) : (A - X = B). (3.54) 


Note that not all B € B can be calculated using any X € Jņ; and all available A € A. The 
set y is not necessarily connected and not necessarily constrained by borders parallel to 
the coordinate axes. Using an enclosing interval X 2 Jņ; will likely create spurious solu- 
tions. 

The controllable solution set applies the same principle to the measurement vector B. 


Definition 3.3 (Controllable Solution Set 5 "4. 


The elements of the controllable solution set are feasible regardless of the chosen point real 
measurement vector B € B. This means there is a suitable regressor matrix A € A for all 
possible point real measurement vectors: 


Y (A B) = (X € R”| (GA € A) A^ (YB € B):(A- X =B)}. (3.55) 


Note that not all A € A can match an element B € B by multiplication with any X € Xy. 
The set 5 ^4, is not necessarily connected and not necessarily constrained by borders paral- 
lel to the coordinate axes. Using an enclosing interval X 2 $y will likely create spurious 
solutions. 
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A very strict criterion is given by the strong solution set. 


Definition 3.4 (Strong Solution set >°\,,) 


The strong solution set includes only parameters X that solve the interval type linear equa- 
tion system for any regressor matrix and any measurement vector: 


Xw (4, B) = (X € R”| (YA € A) ^ (YB € B):(A-X = B)). (3.56) 


None of the limitations of the previous solution sets is necessary for the strong solution. 
Nevertheless, the set $` is not necessarily connected and not necessarily constrained by 
borders parallel to the coordinate axes. Using an enclosing interval X 2 5 w will likely 
create spurious solutions. 

The algebraic solution differs in its definition as it is not quantor based. 


Definition 3.5 (Algebraic Solution Set 5°) 


The algebraic solution is defined by the interval type vectors Xa that solve the interval type 
linear equation system straight forward: 


Y. (A, B) := (X, € IR”| (A - X, = B)). (3.57) 


Even though the elements of the algebraic solution are constrained parallel to the coordinate 
axes, the solution is ambiguous, i.e. there might be several or none solution vectors X, that 
fulfill the equation [Kup95]. 

The different solution sets are related as they are subsets of each other. The united solution 
set is a superset of the algebraic solution set [Kup95], as well as of the tolerable and the 
controllable solution set [Sha96] 


254 (A.B) € 53 (A, B) (3.58) 
Diva (A, B) € D255 (A, B) (3.59) 
Vay (A. B) € $4 (A.B). (3.60) 


The strong solution set on the other hand is a subset of the tolerable as well as of the con- 
trollable solution set [Fie06] 


Xw (4, B) € $4 (A, B) (3.61) 
Dw (A, B) € 3a (A, B). (3.62) 


A visualization of the solution sets and their relations is given in Example 3.6. 
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Example 3.6: 
Consider the following 2 x 2 interval type linear equation system (ILES) taken from [Sha96]: 


(E wa) = (Coa): Bes 
The controllable solution set 5 ^4, and the strong solution set X` are empty for the ILES 
(3.63). It can be seen in Fig. 3.9 that the algebraic solution $^, is a subset of the tolerable 
solution set X ^4 which is a subset of the united solution 5 ^44. It is also clearly visible that 


neither the tolerable nor the united solution can be included in classical intervals parallel to 
the axes without creating spurious solutions. 


Vay =O andYw=0 —4- 


Figure 3.9: Different solution sets for the interval type linear equation system 


The calculation of all given solution sets is computational expensive, as the calculation of the 
hulls is N P-Hard according to [Hor13]. Even to check whether a solution set is empty is still 
an N P-Complete problem as shown by [Sha96]. 


The problem becomes more tractable, if it is regarded from a different point of view. Assume 
a given point real solution candidate X;. The question is now to determine whether the 
solution candidate X, belongs to any of the defined solution sets without calculating the sets 
explicitly. The approach used in this thesis was introduced by [Bee72] and uses the so-called 
theorem of Prager-Oettli [Oet64]. The resulting criterion for interval arithmetic problems in 
Def. 3.6 is used to determine whether X, is a member of the united solution set 5 ^44. 
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Definition 3.6 (Theorem of Prager-Oettli) 
A given solution candidate X, is part of the united solution set 9 ^ 44 i.e. 


Xs € 335 (3.64) 
if and only if X , fulfills the inequality 
de e B.| < A^ |X| + Ba (3.65) 


based on center and radius of the regressor matrix A = (Ac, AA) and the measurement 

vector B = (B., Ba), given by the T interval type measurement values (te) 4 and 
T 

(yx) i1 [Bee72, p. 235]. 


This theorem was extended by [Hla14] to the general A E-solution: 
X, € Dap € |AcX, — Be| € (AR — AX) |X| + BÀ — BY. (3.66) 


A solution candidate vector X, is part of the AE-solution if and only if (3.66) holds. This 
criterion can be specialized to fit the four other solutions sets as given in Tab. 3.1. 


Table 3.1: Conditions for the membership of X to a specific solution set. 


Solution set Condition 

gs (united) [A.X, — B. € AalXs|+ Ba 
doy; (tolerable) |A.X, — B.| € —Aa |X,|+ Ba 
Xay (controllable) |AcXs — Be| < Aa|Xs|— BA 
Dw (strong) |AcXs = B.| < — AA Xs] = Ba 


Further considerations regarding existence and uniqueness of the solution are only avail- 
able for the algebraic solution set »7,. Two approaches for this purpose are sketched in 
Appendix D. 


4 Guaranteed Verification of Point Real 
Systems 


The theoretical foundation of the thesis is developed and illustrated in this chapter. Therefore 
a very simple and comprehensive linear time invariant model structure is used to focus on 
the method itself. The general principles introduced in this chapter can be extended to other 
types of system models. 


4.1 System Setup 


Definition 4.1 (Linear Time Invariant System) 


A discrete time, linear time invariant (LTI) system can be modeled as 


Na Ne 


Uk = ` QiYk—i + 5 CiUk—i (4.1) 
ic j=l 


with the discrete time input u, and output yx, the input and output order Nna and ne as well 
: T T 
as the input parameters |a, a2,...,a4,] and the output parameters [cı,c2,...,Cn.] 


This modeling approach is also known as AutoRegressive system with eXogenous input 
(ARX). 


Based on the model assumption of Def. 4.1 it is possible to set up the specification of the 
nominal system, as given in Def. 4.2. The set of nominal parameters as introduced in Sec. 2.1.1 
is assumed to be given in the specification.* Two possibilities to determine these parameters 
in practice are introduced in Appendix E. 

Throughout this thesis the superscript will be used to denote values that are part of a 
specification or the nominal value. Note that the set of parameters is given by a distinguished 
point real vector for the current LTI setting. 


* 


For other applications, e.g. fault-tolerant control, the method works similarly but the specification is given in a 
different manner. 

The used model assumption does not allow a direct feedthrough as this property is not regarded in the given 
setting. To allow a direct feedthrough the second sum needs to be changed to start from zero, leading to i € 
(0,1,..., nc]. 
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Definition 4.2 (Specification of a Linear Time Invariant System) 


The (direct) specification S} of an LTI system according to Def. 4.1 is given by the mandatory 
values 


e nz, the nominal output order 


e nz, the nominal input order 
T 
. 0" = lat. a3,- Qs Cp Coi i Cuz | ‚the nominal parameters 


and the optional values 


max(n*, 


e You = (Ye) pet me) the initial output values 
e Už a= (uj ana, the initial input values 


leading to the overall specification 


S — 19^, fas es Unit: Yinit } t (4.2) 


If the initial values are known, the future evolution of the system output trajectory is only 
dependent on the input signal. It is thus possible to compare the behavior of the trajectory 
for different inputs. If there are no initial values, the behavior of the trajectory will differ for 
the same inputs in the case of different used initial values. If they are provided, there need 
to be at least kmin = max (nž,nž) + 1 initial values to enable the first evaluation of the 
autoregressive system description according to Def. 4.1. 


It is assumed that the nominal system is developed and built and ready to be verified. Thereby 
the verification object (VO) is assumed to be available as a (physical) black box with one or 
more input and output ports. It is possible to excite the system via the input and to measure 
the resulting output. Further insights, like internal structure, components and wiring, soft- 
ware, plans or internal states are not accessible. This approach can be applied in various states 
of system development. Therefore there is a wide range of exact physical representations of 
the VO black box such as models, program code, components or units. 


The verification method is running on a digital device that not necessarily generates the input 
signal itself. Therefore input and output values need to be measured to be available for the 
verification method. Measurement data is always subject to noise which is assumed to be 
modeled throughout the thesis based on the following definition. 
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Definition 4.3 (Sensor Noise Properties) 

All available information about the VO is given in terms of measurement data ofthe input 
N and output TEUER. ae The measurement data is obtained using sensors 
providing guaranteed notions of sensor precision that allow interval type enclosure of the 
measurement values. 

All measurement values Umeas,k and Ymeas,k are extended to intervals such that the true 
system values Utrue,k and Ytrue,k are guaranteed to be included in the interval Umeas,k 
and Ymeas,k, respectively. There are three different ways these guarantees can be obtained: 


a) Absolute deviation 
The sensor precision is denoted by a maximum deviation of +6“. This leads to the 
interval enclosure of 


Utrue,k € Uk = [mens ie m ns meas ,k + Onl (4.3) 


Utrue,k € Yk = | esed E Bi Umeas,k + ôs] : (4.4) 


b) Relative deviation 
In this case the sensor precision is given in terms of a relative deviation of ô” € [0, 1] 
which is leading to 


Utrue,k € Uk = [Umeas,k Jj (1 = 1); Umeas,k * (1 "E or) (4.5) 


Utrue,k € Yk = ae: . (1 = Ôr) Ymeas,k ` (1 EE 87)] . (4.6) 


c) Combined deviation 
A common case is the combination of the both aforementioned deviation types. The 
deviation is defined to be 6” € [0, 1] times the current measurement value but at 
least £0", resulting in 


: ' ĝe 
Utrue,k € Uk = Umeas,k * [min (a — 01), (1 = 9) 
Umeas,k 


max (a +8"), (1 i x -)) (4.7) 

Umeas,k 

ó2 
Umeas,k ' [min (a = M (1 = 23) : 
Ymeas,k 
6° 

max (urn, (+) , (4.8) 

Ymeas,k 


| 


Utrue,k € Yk 


The properties of Def. 4.3 are used to set up interval type enclosures of the measurement data 
that are guaranteed to include the true system value. The resulting structural setup of the 
measurement process is depicted in Fig. 4.1. 
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Figure 4.1: Structure of measurement setup 


A basic assumption in the field of parameter identification is the property of persistent ex- 
citation according to Assumption 2.1. The information provided in any data set is highly 
dependent on the input signal that was used to generate the output. According to [Ast95, 
p. 63ff] there are several methods to ensure persistent excitation of a system. Exemplary per- 
sistently exciting inputs are e.g. white noise, pseudorandom binary sequences or a moving 
average process [Ise10, p. 251]. One possible excitation procedure fitted to the specific settings 
regarded in this thesis was developed in [Rie17]. The main idea is sketched in Appendix F. 


4.2 Time Invariant Full Consistency 


The general notion of consistency introduced in Chapter 2 is now transferred to the specific 
setting of LTI systems. The direct specification $7 includes one distinctive point real nominal 
parameter vector. Thus all results show F'ull Consistency according to Def. 2.2. 

The observed behavior of the regarded VO is given in terms of input output measurement 
data. The nominal behavior is specified according to Def. 4.2. The VO is called full consistent 
with its specification if the measurement data can be explained by all specified parameters. 
The verification question is formally stated in Problem 4.1. 
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Problem 4.1 (Time Invariant Full Consistency) 


Is the nominal system, specified by a direct specification 
Sa = 18"*,nj er in}; (4.9) 


full consistent with the input-output behavior given by the interval type enclosures of T 
measurement values 


[U ena; Y meas] = eee Cen P (4.10) 


i.e. can the measurement data be explained by the nominal system? 


Problem 4.1 can be solved using the united solution set according to Def. 3.1. 


Proposition 4.1 (Time Invariant Full Consistency) 


The interval enclosure of the measurement data [Umeas, Ymeas], given for the dis- 
crete sampling points k = {1,2,...,T}, leading to the interval type regressor ma- 
trix Ameas and the interval type measurement vector Bmeas, of a VO is called 


full consistent with a direct specification Sž, if the specified system parameters 
T 
o= lat. 05,..., Qs, C], C... Cn. | are part of the united solution set gs given 


by the measurement data, i.e. 


(O* € 3055 (Ameas; Bmeas)) € Full Consistency” => Full Consistency. (4.11) 


Proof: 
gi 


The nominal parameter vector O* — lai. 03; «One i OL Peg A , given by a direct 
specification S} according to Def. 4.2, can be interpreted as a solution candidate for the ILES 
(3.46) which is set up by the interval type enclosure of the measurement data. If O* is part 
of the solution set of the ILES (3.46), the specification 57 is able to explain the measurement 
data. 

The problem is formulated in Kaucher interval arithmetic, therefore it is necessary to define 
which solution set is used. Considering the interval enclosure of the measurement data given 
in Def. 4.3, it is obvious that it is not possible to determine the true value Utrue,k and Ytrue,k 
as there are two distortion steps between the true values and the interval enclosure. First 
the true value is changed by the measurement random noise e. Second, the sensor is only as 
precise as given by its property. 
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However, it is guaranteed that the true values of Utrue,k and Ytrue,k are included in the mea- 
surement intervals 


Utrue,k € Umeas,k (4.12) 


Utrue,k € UYmeas,k- (4.13) 


Starting from time kmin = max (nz, nz) + 1 there are enough measurement values to set 
up the regressor equations. Each additional measurement value leads to an additional row in 
the regressor matrix A. 

The unknown true values can be assumed to form an unknown true point real regressor 
matrix 


Yrue,kmin—l ze Yırue,kmin—nı Utrue,kmin—1 QS Utrue,kmin—n* 
Ytrue,kmin ud Ytrue,kmin+1—n* Utrue,kmin a Utrue,kmin+t1—n* 
Atrue = 
Utrue,T—1 "Ut Ytrue,T—1—n* Utrue,T—1 Ut Utrue,T—1—n* 
a c 
(4.14) 


and an unknown true point real measurement vector 


Birue = forme Ytrue,kmin+1» -> „Ytrue T] . (4.15) 
These elements are linked via an unknown true parameter vector Otrue, that fulfills 
AtrueQtrue = Birue- (4.16) 
Based on the enclosure of the true values in Def. 4.3 holds: 


Atrue € Ameas (4.17) 
Birue € Beas. (4.18) 


With the set definition (3.53) follows that ©;,,- is an element of the united solution set 
Yaa (9A sens Pacis): 

It is impossible to determine the true values A;,,, and Birue from the given measurement 
data. Therefore each point real element of the interval type regressor matrix and the interval 
type measurement vector is a possible true value. Thus the whole united solution set can be 
considered as correct solution of the ILES. 

The given direct specification S} shows Full Consistency” with the given measurement 
data [Umeas; Ymeas|, if and only if the nominal parameter vector O* is part of the united 
solution set 5 ^44 (Amcas; Bmeas). Due to the underapproximating property holds 


Full Consistency” > Full Consistency (4.19) 


and thus the VO is full consistent in the sense of this thesis. 
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The calculation of the united solution set is computationally expensive as introduced in Sec- 
tion 3.1.2. However, it is not necessary to calculate the whole solution set in this setting as 
there is a candidate solution given in form of the specification. Thus it is sufficient to check 
whether the specified parameter vector ©* is part of the united solution set without calculat- 
ing the solution set explicitly. Prop. 4.1 can be checked very efficiently using the theorem of 
Prager-Oettli according to Definition 3.6 by evaluating the single equation (3.65). Therefore 
(3.65) is reformulated for the given measurement values to 


|Ameas,cO” = Bmeas,c| < Ameas,A |O*| + Bmeas,A e0*c Mss (Ameas: Bmeas) . 
(4.20) 


As stated in Section 3.1.3 the existence and uniqueness of the solution sets is still an open 
question. A necessary condition for the existence of any solution set is that the ILES (3.46) is 
solvable. For the general overdetermined setting, this can be checked using the approaches 
given in Appendix C. However, their application is limited as the problem is in general N P- 
hard. The only further considerations regard the algebraic solution set and are sketched in 
Appendix D. 


Note that the introduced method preserves time-invariance when checking for consistency. 
This is due to the used specification and represents a main difference to the direct image based 
methods used in fault detection as introduced in Section 2.2. This property will become even 
more clear in Chapter 5. 


A further necessary condition is persistent excitation of the VO which is given in Assump- 
tion 2.1, developed to ensure full rank according to Assumption 3.1. 

Therefore it is in general not guaranteed that there is a nonempty united solution set avail- 
able and thus there are situation in which the proposed method is not applicable. 

However, it is possible to facilitate a favorable situation by proper experiment design. One 
possibility to determine a beneficial excitation signal is given in Appendix F. 


The application of time invariant full consistency for LTI systems is demonstrated in Exam- 
ple 4.1 and was presented to the scientific community in [Sch17b] and [Sch17c][Sch19]. 
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Example 4.1: 


This example shows the verification of a linear, time invariant system as introduced in 
Prop. 4.1. Assume the following direct specification 

Si1 = 19* = [0.9825, 0.0675] ,n, =1,n2 = 1, Unit = [0], Vine = 0}: (421) 
The simulations are done with a virtual VO, correctly implemented as discrete time linear 
ARX system 


Yr = 0.9825yy 1 + 0.0675u,—1 (4.22) 


with sampling time At = 1s. The system is excited using a noise signal with uniformly 
distributed amplitude Utrue,n € [0, 10] with mean Umean = 5. It is assumed that the input 
is measured using a sensor with a maximum relative fault of 0j, = 0.05. The resulting 
enclosure of the input measurement signal Umeas is depicted in the first subplot of Fig. 4.2. 
Nevertheless the system (4.22) is fed with the undisturbed input signal Utrue. The resulting 
output signal is measured using a sensor with the same properties as the input sensor. The 
enclosed measurement output signal Y meas is depicted in the second subplot of Fig. 4.2. 
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Figure 4.2: Measured input signal Umeas with 67, = 0.05 
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It can be shown that Prop. 4.1 holds for $5, and the measurement data Yeas and Umeas 
which proofs full consistency of measurement and specification formally. 

A graphical representation is given in Fig. 4.3. The lines depict borders of the united solution 
set, generated by the different rows of the regressor matrix. Feasible parameters need to be 
located in between the borders of all rows of the measurement matrix. The parameters given 
in specification 97 , are marked with a green cross and form a feasible solution of the given 
problem as they are located within the united solution set of all measurement data. Hence it 
is possible to explain the measured data with the parameters given in specification S% ,. 

All given measurement values are used in this example to set up the regressor matrix and the 
measurement vector. This leads to Ameas € IIR9*2, In this case, the full rank Assump- 


tion 3.1 for n = 2 parameters leads to rank (Ameas) = 2 which is fulfilled for the given 
dynamic and excitation signal. 

An example with failed verification can be given for the case that the verification method is 
applied using a different specification on the same measurement data. For this purpose 


Sao (0* = [1.15, 0.08] , n? = 1,n; = 1, U; [0] ,Y; (0]) . (4.23) 


init — >t init Z 
is used. 
A graphical representation of the parameters is given by the red mark in Fig. 4.3. In this 


case, the parameters do not lie within the borders of the united solution set of the given 
measurement data. Thus the system is not verified. 
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Figure 4.3: Visualization of the united solution given by the measurement data 
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4.3 Conclusion 


The main idea of Kaucher arithmetic based verification was introduced in this chapter. There- 
fore the precise system and problem setting was defined and explained. The problem setup 
leads to a very distinct situation with only poor knowledge about the true measurement val- 
ues which are enclosed in intervals. If this setting is regarded from a different point of view, 
it can be interpreted as a specific quantor based solution set definition that matches exactly 
the united solution set introduced in the mathematical preliminaries. This property can be 
used to check full consistency of the measurement data and the specification in a set member- 
ship procedure that is computationally very effective. The assumption of a full rank interval 
regressor matrix leads to preliminaries on the sensors and the noise assumptions. It is pos- 
sible to check whether a specific solution candidate is part of the united solution given by 
measurement data. This property was demonstrated using an illustrative example. 


The main advantage of the introduced method is that it focuses on the united solution set. 
Thereby it is possible to avoid wrapping and dependency effects and to calculate a solution 
set free of spurious solutions. This property is very beneficial in the case of safety critical 
systems as it avoids type II errors (hidden alarms). 


5 Guaranteed Verification of Interval Type 
Systems 


The basic idea introduced in the previous chapter is now extended to an interval type speci- 
fication. Thus the parametrization is given by an interval type vector O* instead of a point 
real vector O*. 


Definition 5.1 (Interval Type Specification of a Linear System) 


An interval type specification S7 of a linear system is given by the mandatory values 
e nz, the nominal output order 


e nz, the nominal input order 


. ©* = lai. 05, ..., 05,4, Ci C3, -+ Cp» |, the interval type nominal system param- 
eter vector 


and the optional values 


-Y* = ( ee the initial output values 
e Usa = (uj) maximam) the initial input values. 


This leads to the overall specification 


oj = 1O* , as Uim Yoal- (5.1) 


Ti 


Based on this specification the system is implemented. It is assumed that the resulting VO is 
given in a form that provides the input and output signals as described in the specification. 
Again, this can be the case for a variety of test objects, depending on the specific point of 
the development cycle for which the specification 57 was defined. Methods to determine the 
nominal parameters in practice are given in Appendix G. 


Even though the specification is now given by interval type values, the implemented system 
has to provide real output data at any given time. Thus the real implementation of the VO 
has to use a specific real parametrization. This leads to the definition of interval type linear 
systems as given in Def. 5.2. 
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Definition 5.2 (Interval Type Linear System) 


A discrete time, linear, interval type system can be modeled as 


Na Ne 
Vk = 5 Qi kUk—i + 5 Ci kUk—i (5.2) 
i=l i=l 


with the discrete time input u, and output yx, the input and output order n, and n, as well 
as the time variant parameters 


I: 
Op = [a1 k; 02,5. n, jo CL Co nah € ©. (5.3) 


The necessary data is available as disturbed, discrete time measurement data enclosed in 
intervals according to Def. 4.3. Also the persistent excitation Assumption 2.1 and the full 
rank Assumptions 3.1 are still required to hold. 


The given system definition leads to the time variant regressor vector 
Ak = [Yk—15 Yk—-23 +++) Ukenas Uk—1; Uk—2; ++) kn] (5.4) 
and thus the system equation can be transferred to 
Uk = AkOk (5.5) 


for a specific time step k > kmin with kmin = max (na, nc) +1. 


This can be interpreted as the realization of a time variant system whose interval type spec- 
ification is given according to Def. 5.3. 


Definition 5.3 (Interval Enclosure of Time Variant Parameter) 


The parameter values ©; evolve during a specific time k € {1,2,...,T} and can be en- 
closed in the interval 
T 
Oc O= a, a), " na] (5.6) 


i GT (NT 
with n. = na + n, and OU) = [nin (C ). ) ,max & ^ ). denoting the 
;—1 ;—1 


minimum and maximum value of the i-th component within the regarded time. 


The time variance is given only in the parameters, the model structure, especially the model 
orders n, and ne, are time constant. Furthermore this interpretation is not necessarily ben- 
eficial for all time variant systems as the resulting interval enclosures can be very large de- 
pending on the time variant dynamic of the system parameters. 


5.1 Interval Type Full Consistency 


5.1 Interval Type Full Consistency 


In the case of an interval type specification both consistency definitions (Full Consistency 
and Basic Consistency) according to Def. 2.2 and Def. 2.1 are possible. In this section, the 
idea of Full Consistency is extended to interval type systems. Then the situation is relaxed 
to Basic Consistency in the next section. 


Problem 5.1 (Interval Type Full Consistency) 
Is the nominal system, specified by an interval type specification 


S; = 1e*, Tis fic, Units Yol ? (5.7) 


full consistent with the input-output behavior given by the interval type enclosures of T 
measurement values 


[U mens Vincas| = OR AEEA = (5.8) 


i.e. do all elements of the parameter vector ©* fulfill Prop. 4.1 for all measurement data? 


The full consistency problem is solved by Prop. 5.1. 


Proposition 5.1 (Interval Type Full Consistency) 

The interval enclosure of the measurement data [Umeas, Ymeas| given for the discrete 
sampling points k = {1,2,..., T}, forming the regressor matrix Ameas and the measure- 
ment vector Bmeas of a VO, is called to be full consistent with an interval type specification 
S1, if the complete set of specified parameters ©* is part of the united solution set X`- 
given by the measurement data, i.e. 


(O* C Sas (Ameas; Bmeas)) © Full Consistency” => Full Consistency (5.9) 


Proof: 
Full consistency follows directly from applying Prop. 4.1 to all possible point real parameter 
vectors O* € ©* given in the interval type specification $7. 
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The inverse relation of Prop. 5.1 leads to the implications given in Prop. 5.2. 


Proposition 5.2 (Inverse of Full Consistency) 


If there is at least one O* € ©* that does not show full consistency according to Prop. 4.1, 
the interval type specification S7 is not full consistent with the measurement data. 


Proof: 
Full consistency according to Prop. 5.1 is defined for all parameters O* € ©*. A parameter 
O* € ©* that does not show full consistency according to Prop. 4.1 leads to 


iu g y» (Ameas: Bineas) : (5.10) 


and thus the interval type specification 57 is not full consistent with the measurement data. 


Theoretically, Prop. 5.1 can be checked by applying the verification equation (4.20) based on 
the theorem of Prager-Oettli to each parameter vector O* € ©*. However, to realize this 
approach in an algorithmic implementation it is necessary to draw nNeneck discrete samples 
from the continuous intervals. The number of discrete parameter vectors to check neheck 
thus becomes a relevant design parameter. According to Prop. 5.2, a single inconsistent vec- 
tor falsifies the full consistency property. This leads to the requirement that neheck has to be 
very large to cover the given parameter range sufficiently. Even though a single evaluation 
of the theorem of Prager-Oettli is computationally very effective as stated in Chapter 4, the 
computation time rises proportionally with ncnecr- Additionally, as this thesis aims on cal- 
culating guaranteed results, the step size used to sample the interval type parameter vector 
needs to be very fine, even tending to zero. This high resolution needs to be applied to each 
component of the parameter vector. Afterwards it is used to build all possible combinations 
including the samples of the different components. Assuming a resolution of n, samples on 
each component of O* leads to 

check = ng" ne (5.11) 
applications of the theorem of Prager-Oettli . This number increases polynomial with n, and 
leads to large computation times for sufficiently high resolutions. Thus the sampling based 
approach is computationally infeasible. 


Restructuring the problem can be used to avoid the necessity to cover the whole parameter 
area. The verdicts can be calculated based on the vertexes Y ofthe nominal parameter set only 
and then can be generalized to the whole nominal set if convexity properties are fulfilled. 
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The maximum number of points to check is thus reduced to 
Ncheck = niin: (5.12) 


which resembles a computationally feasible number, especially for low system orders n; and 
nz. The vertexes V ofthe hyperrectangle given by the interval type parameter vector O* can 
be determined according to Def. 5.4. 


Definition 5.4 (Vertexes of a Hyperrectangle) 


The nominal interval vector @* € IR"«*"** defines 

n = Patni (5.13) 
vertexes V € IR" +"! of a hyperrectangle that can be indexed using a decimal index 
Vaec € (0,1,..., m — 1}. The index is subsequently transformed to its binary representa- 
tion Vy; that can be interpreted as (1 x n7 + n) dimensional vector were the i-th vector 
component is denoted as ve. 
The specific values of the vertexes V,,,.. can be generated by interpreting the binary in- 
dex Vpin component wise for i € {1,2,...,n% + nz) and extracting the limits from the 
respective nominal parameter vector element: 


*(i) @) _ 
Den Um (5.14) 
“lem gines ed 
An illustration of Def. 5.4 is given in Example 5.1. 
Example 5.1: 
Consider the following (2 x 1) interval vector with n, = nz = 1 
e* = [[2, 3], [4, 6]” (5.15) 


The resulting rectangle has n = 2? = 4 vertexes with vaec € {0,1,2,3} and i € {1,2} 
leading to the indexes given in Tab. 5.1. 


Table 5.1: Vertexes of a hyperrectangle 
Decimal index Vgec | Binary index Vy; | Coordinates V,,,... according to (5.14) 


0 


[2 4] 
[2 6] 
[3 4] 
[3 6] 


w Ne 
S) 
IAS 
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It is now possible to set up an alternative formulation of Prop. 5.1, that solves Problem 5.1 
with vertex based full consistency. 


Proposition 5.3 (Vertex Based Full Consistency) 

The interval enclosure of the measurement data [|Umeas, Ymeas| given for the discrete 
sampling points k = {1,2,...,T}, forming the regressor matrix Ameas and the measure- 
ment vector Bmeas of a VO, is called to be full consistent with an interval type specification 
S*, if all vertexes V defined by the set of specified parameters O* are located in the same 
orthant and are part of the united solution set 544 given by the measurement data, i.e. 


(Y € ^34 (Ameas; Bmeas)) & Full Consistency” > Full Consistency. (5.16) 


Proof: 
The united solution set can form various shapes, but it was shown by [Sha10] that the general 
AE-solution is convex within each orthant. As the united solution set is a special case of the 
general AE-solution, this property does also hold for 755 (Ameas; Bmeas): 

The specified parameters ©* and thus the resulting vertexes V are all located within the same 
orthant. 

The theorem of Prager-Oettli can be checked for the n = 2"2+"< vertexes in finite time. 
Based on the direct application of the definition of a convex set given in [Bro08, p. 662] 
follows: 

If (3.65) holds for any two of the vertexes V; and V}, with i, j € {1,2,...,n — 1}, i Æ j, all 
vectors © = AV; + (1 — A)V;, with 0 € A < 1 are also part of the united solution set. 
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Example 5.2: 


Assume the same setting as in Example 4.1. However the specification is now given as an 
interval type specification 


S= [e = [[0.9725, 0.9925], [0.0665, 0.0685]]" , nz = 1, nt = i (5.17) 


The simulations are done using the linear discrete time ARX system (4.22) leading to the 
same measurement data as given in Fig. 4.2 in Example 4.1. The resulting verification setting 
is depicted in Fig. 5.1. The nominal parameters given in S7 are depicted as green square. It can 
be seen that all four vertexes V = (Vo, Vi, V2, V3} are located within the united solution set 
given by the measurement data. The specification and the measurement are thus guaranteed 
to be full consistent according to Prop. 5.3. 
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Figure 5.1: Example setting for an interval type specification S* 
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5.2 Interval Type Basic Consistency 


Until now the set of VO behavior was assumed to be considerably larger than the specifica- 
tion. When using interval type specifications this is not necessarily the case. It is possible 
that the specification set is of the same size as the set of VO behavior or even larger. Therefore 
it is not longer possible to enclose the whole specification in the VO behavior. The resulting 
verification question can be formulated as follows: 


Problem 5.2 (Interval Type Basic Consistency) 
Is the nominal system, specified by an interval type specification 


5; = 197 nas Ne, U; it: Y , (5.18) 


Ti 


basic consistent with the input-output behavior given by the interval type enclosures of T 
measurement values 


[Us eas; Y ineas] = ee 13 CE (5.19) 


Le. is there at least one parameter vector O* € ©* that fulfills Prop. 4.1? 


The Venn chart of a basic consistent setting is depicted in Fig. 5.2. Due to the inner enclosure 
of the VO behavior the achieved verdict is still type II error free. The consistent set is given 
by the green shaded square. 

This set it formally stated in Prop. 5.4 and solves Problem 5.2. 


Nominal Behavior 


Consistent Set 


VO Behavior IA 
chavo 


Underapprox. 
VO Behavior 


Figure 5.2: Basic consistent result for a large specification 
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Proposition 5.4 (Interval Type Basic Consistency) 


The interval enclosure of the measurement data [Umeas, Ymeas| given for the discrete 
sampling points k € (1,2,..., T), forming the regressor matrix Ameas and the measure- 
ment vector Bmeas of a VO, is called basic consistent with an interval type specification 
S*, if there is a nonempty consistent set, i.e. a nonempty intersection between the nominal 
set O* and the united solution set 9 ^44 (Ameas; Bmeas) 


(O*A Y a4 (Aneas; Bmeas) # 0) €& Basic Consistency” > Basic Consistency. 
(5.20) 


Proof: 
All parameter vectors included in the interval specification ©* are suitable representations 
of the correct system behavior. Thus 57 can be interpreted as a set of direct specifications 
S7. Each of these direct specifications 57 € 57 can be used to check consistency according 
to Prop. 4.1. If there is at least one full consistent direct specification included in the interval 
specification, the VO behavior can be explained by this parameter and the VO is denoted as 
basic consistent. 


Basic Consistency” is sufficient for the genuine system, as that there is at least one param- 
eter within the nominal set that is able to explain the measurement data. 


5.2.1 Algorithmic Solutions 


A straight forward approach to check basic consistency uses the vertexes only. However, 
the shape of the resulting consistent set will change if one or more vertexes are inconsistent. 
This change is depicted exemplary for the 2D case in Fig. 5.3. 


x x x x 


x x x 
Feasible 
Points 


Figure 5.3: Degradation of the consistent set for different consistent vertexes (green) 


The shape changes from a rectangle, in case all four vertexes are part of the consistent set, 
to a single point if only one vertex shows consistency. The main drawback of this procedure 
is that there might be a consistent set, even if no initial vertex is consistent. This situation is 
exemplary depicted in Fig. 5.5. 
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Also, the basic shape of the resulting consistent set changes which can be a disadvantage for 
the further algorithmic processing of the result. 

A solution for the first problem is given by checking more points that are not a vertex. How- 
ever this leads to the sampling based approach as introduced in the previous section, with 
the respective runtime limitations explained there. 

The problem of finding the right resolution in the sampling based approach can be solved 
using optimization methods. The idea is to use the logic of optimization algorithms to guide 
the sampling process.° Optimization procedures can be used to determine which discrete 
points to check if one or more vertexes of the specification are not element of the consistent 
set. 

To use optimization methods for the choice of sampling points the problem can be refor- 
mulated to a feasibility problem. The interval type measurement data is used to set up the 
constraints that frame the united solution set 5 ^ 44. As derived in the previous chapter, every 
parameter within the united solution set 5 ^44 is able to explain the measurement data. The 
constraints limit the search area of the feasibility problem. Feasible solutions need to fulfill 
all constraints. Parameter values that are located outside the united solution set 5 744 are not 
feasible with respect to the constraints. 


The interval type specification 57 includes the nominal parameter vector ©* which can also 
be denoted as nominal set M. The nominal set represents the maximum area of potentially 
consistent parameters. This initial restrictions can be stated in terms of linear inequality 
constraints as follows: 


Definition 5.5 (Constraints Given by the Nominal Set) 


The interval type parameter vector ©* given in the specification S7 can be used to set up 
the set cy of 2(n% + ni) linear inequality constraints that restrict the feasibility problem 
to the nominal set N: 


c (e) = GY 9) <9 (5.21) 
creme) (0) m 7) d: g? < 0 (5.22) 
with the number of parameters i € {1,2,...,n*% +n*}. 


Suitable optimization methods are e.g. grid search, golden section search or dichotomous search [Wil64] if there 
is only minimal information available. If there is further knowledge about the shape of the problem, there are 
more sophisticated algorithms that can direct the search effort very efficiently into the relevant regions, e.g. 
Newton method, simplex method or interior point method [Noc06]. 
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The united solution set 5 744 can be reformulated in terms of the interval type measurement 
data as linear matrix inequality constraints (LMI) as given in Def. 5.6. 


Definition 5.6 (Constraints Given by the Measurement Data) 


P T : 
The measurement data [U cas; Y eas] = [rs m , (Ymeas,k)p—1| forming the 


regressor matrix Amcas and the measurement vector Bmeas of a VO, can be used to set up 
the set cm of 2(T — max (n*,n*)) linear inequalities 


VIC) = —(Ameas : 9) + Bil... <0 2) 
a ie (6) = (Ameas ` oj” = E < 0 (5.24) 
with the number of rows in the regressor matrixi € (1,2,..., T' — max (nž,nž)}, ie. the 


number of system equations instantiated by different measurement points. 


Each pair of constraints represents the upper and lower bound of the solution set derived 
from one line of the ILES (3.46). They can be interpreted as hyperstripes in the parameter 
space, leading to the setting shown in Fig. 5.4 for the 2D case. Using constraints according to 
Def. 5.6 limits the search area of the optimization algorithm to the inner approximation and 
thus guarantees that there are no type II errors possible in the resulting consistent set. The 
consistent set is given by the shaded region. 


Figure 5.4: Exemplary constraints of the feasibility problem in 2D 


Note that the number of constraints in c 1 is directly related with the number of used sam- 
pling points T and that the number of parameters n = m; + n; is of minor importance. 
However there have to be at least T = kmin = max(nž,nž) + 1 samples in the measure- 
ment to set up the first hyperstripe. It is now possible to determine an alternative solution 
for Problem 5.2, based on the constraints of Def. 5.5 and Def. 5.6. 
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Proposition 5.5 (Feasibility Based Basic Consistency) 


The interval enclosure of the measurement data [|Umeas, Ymeas| given for the discrete 
sampling points k € (1,2,..., T), forming the regressor matrix Ameas and the measure- 
ment vector Bmeas of a VO, is called basic consistent with an interval type specification 
S; , if the consistent set C is nonempty, i.e. if there is at least one solution O that fulfills all 
constraints 


((ev(6) < 0) ^ (e(6) < 0)) < Basic Consistency” = Basic Consistency 
(5.25) 


Proof: 

According to Prop. 5.4 there needs to be at least one parameter that is part of the united so- 
lution set 5 ^44 as well as part of the parameter set given in the interval type specification 
S; to ensure basic consistency. All parameters that fulfill the constraint set cm of Def. 5.6 
are part of the united solution set 5 ;;. All parameters that fulfill the constraint set cy of 
Def. 5.5 are part of the parameter set given in the interval type specification 57. If there is at 
least one parameter vector © that fulfills cy and cm, all conditions for basic consistency are 


fulfilled. 


The resulting situation is exemplary depicted in Fig. 5.5. It can be seen that no initial vertex 
shows consistency, as none of them is part of the underapproximation. Therefore there is no 
vertex based basic consistency. Feasibility based basic consistency can be achieved as there 
is a consistent set within the inner approximation and the nominal set. The consistent set is 
depicted as the shaded green area. 
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Figure 5.5: Example setting showing feasibility based consistency 
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5.3 Conclusion 


This chapter introduced the extension of the Kaucher based method to interval type specifi- 
cations. The resulting nominal set is required to be located within one orthant. The definition 
of full consistency is still applicable if all parameters given in the interval type specification 
are part of the united solution set. A straight forward approach was introduced that used the 
orthant wise convexity of the united solution set to define full consistency based only on the 
vertexes of the specified parameter set. 

Nevertheless, full consistency is a rather strict criterion and it is likely that there are some 
specified parameters that are part of the united solution and some that are not. Following 
the concept of basic consistency, it is sufficient to show that there is at least one parameter 
vector that is part of the united solution set as well as of the specified parameter range. This 
property can be verified by checking individual arbitrary points for consistency. The choice 
of these points can be structured by using the concept of a feasibility problem. Therefore the 
notions defining the united solution as well as the specified parameter set are transferred to 
linear matrix inequality constraints. If there is a solution of the feasibility problem, the VO 
and the specification are guaranteed to be basic consistent. 

The method is still free of type II errors as the feasibility based consistent set is constrained 
by the genuine united solution set. 
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Switched hybrid systems consist of two distinct parts with different properties and modeling 
goals. The dynamic part is used to model the plant dynamics as introduced and used in the 
previous chapters. The additional discrete event part models the superimposed switching 
logic. Based on logical rules, the discrete event part can activate different discrete states that 
will show different dynamic behavior. Different operation modes can thus be modeled as 
several subsystems, showing individual behavior. All subsystems are interconnected by the 
discrete event switching mechanism. 

The switched hybrid system structure is depicted in Fig. 6.1 and is formally introduced in this 
section’. 


State Machine Z 


: Switch signal / State Signal 


Set of Dynamic Systems S 


Hybrid System H ı 


a Se eee ee a a a a a nn m 


Figure 6.1: Structure of the hybrid system model H 


The hybrid system H consists of a set of dynamic subsystems S and a superimposed switching 
mechanism represented by the state machine Z. 


The subsystems s(? € S show different behavior based on an individual parametrization. 
The state machine produces a switch signal which resembles its current discrete state. This 
signal is used to control an input and an output switch that determines which subsystem is 
activated. 


7 To improve readability the term “hybrid system” is used instead of “switched hybrid system" throughout this 


thesis. 
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The activated subsystem is fed with the general input signal and the resulting output signal 
is connected to the output of the hybrid system. It is possible to use the optional input and 
output values that are specified in the nominal system to start the active subsystem after a 
switch. Otherwise the current input and output values are kept across the switch. The input 
and output signals are also fed to the state machine where they are used to update the discrete 
state of the discrete event system. Due to this extended structure there are additional subjects 
included in the verification question that will be covered in this chapter. First the formalized 
model structure needs to be appended to a hybrid formulation. Therefore the discrete event 
part is modeled in the following, based on [Cas99, p. 66ff]. 

Then the verification problem is split in two subproblems: verification of the dynamic part 
and verification of the discrete event part. If both parts are verified individually, the next 
step is to examine their connection and interaction. This is done in three steps of increasing 
complexity. 

Initially, the setting is simplified by assuming the discrete state to be measurable. This setting 
is used to introduce the basic hybrid method. Second this assumption is dropped such that 
the current active discrete state needs to be determined for a given set of switches. Finally an 
algorithm is developed that is able to determine the switching times and the active discrete 
states from measured input and output data only. 


The necessary knowledge for the verification procedure is thus the same as in the previous 
chapters except that there are now several nominal systems and an additional specification 
of the discrete event system part. 


Definition 6.1 (Discrete State) 


A discrete state 
q® € Q:= faa, T qo) (6.1) 


is a vertex of a graph e.g. of a state machine. It is part of a given set of discrete states Q. 


Note that for the ease of notation the “discrete state" is called “state” in the reminder of this 
thesis. 
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To implement logical conditions and constraints in the state machine a set of events is de- 


fined: 


Definition 6.2 (Event) 


There is a set of events 
e fen e@,... e0}, (6.2) 
Each event e) is dependent on specific activation limits 
i) = e, r] , (6.3) 


with i € {1,2,...,n.}, defined for a given enabler signal W := (wk)1 The event e? 
is defined to be active as long as the value of the enabler signal W lies within the given 
thresholds 


wy, € IM. (6.4) 


Note that Def. 6.2 does not pose any conditions on the activation limits. Therefore it is pos- 
sible that several events are active at the same time. In the context of this thesis, events are 
allowed to be active for several time steps, e.g. as long as the enabler signal stays within the 
specified limits. 

The transitions of the state machine are defined by a transition function. 


Definition 6.3 (Transition Function) 
A transition function 


f:QxE>Q (6.5) 


represents a directed connection between two states, labeled by an event. In general f is a 
partial function on its domain. 


The notation t® : f (q®, e) = q) means that transition t(!) forms a directed connec- 
tion from state q(U to state q®), dependent on event ei), 

A transition can change the state of a state machine if the assigned event is active, but not 
necessarily has to. This is due to the fact that several events can be active at any given time, 
but there is not more than one transition allowed to conduct a switch. However it is also 
possible that the state does not change even though there are several activated transitions. 
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Based on the above definitions it is possible to set up the state machine representing the 
discrete event part. 


Definition 6.4 (State Machine) 
A state machine is defined to be given by the 4-Tupel 


Z := foe, fa), (6.6) 


with a finite set of states O, a finite set of events £, a transition function f and an initial 
state qh), 


The state of the discrete event part can be used to determine the system orders and parame- 
ters of the dynamic part necessary to set up a system according to Def. 5.2. This leads to the 
definition of a state dependent, discrete time, linear, interval type system: 


Definition 6.5 (State Dependent Discrete Time Linear Interval Type System) 


The state dependent, discrete time, linear, interval type system can be modeled as 


Na (qk) Ne (dk) 
sit) = y, = — 5 ai(dk)Yk—i + A ci(qx)ui—i (6.7) 
i=1 i=1 


with the discrete time input uy, output yy and state qu. The input and out- 
put orders na(qx) and n.(gr) as well as the interval type system parameters 
© = [ai (qk), a2(9r),-- - ;An.(q,) (dk); €1 (dk); 2 (dk). +++ €ne(qy)(Ge)] are dependent 
on the current state qp. All subsystems s(4*) form the set of subsystems 


S= fa 0, stah, (6.8) 


Each dynamic subsystem is directly linked with a discrete state. Therefore the state depen- 


dent dynamic subsystem s(0®) is denoted by s(? for the ease of notation. It is now possible 
to define the overall hybrid system. 
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Definition 6.6 (Hybrid System) 
A hybrid system H consists of two system parts: 


Discrete Event Part The superimposed switching mechanism given by a state machine 


Z - (9,5, fa}, (6.9) 
according to Def. 6.4. 


Dynamic Part The discrete time linear interval type systems are given by a finite set of 
subsystems 


S= fa s, s (6.10) 


where each subsystem s(4*) is active if and only if qj, is the current state. The sub- 
systems s(4*) are defined according to Def. 6.5. 


In general it is not possible to measure the current state of the state machine. However, for 
didactic reasons Assumption 6.1 is introduced and later dropped. 


Assumption 6.1 (Measured State Signal) 


The current state of the state machine can be measured and it is correctly given in the state 
signal 


T 
Qmeas = N G (6.11) 


The measured state signal Qmeas consists of several segments with different active states. A 
change in the active state is called switch. 


Definition 6.7 (Switch) 


The first time index k at which a new state is active i.e. 


dmeas,k £ (Qmeas,k—15 (6.12) 


is called switch k+. 


The switches within a state signal are additionally indexed in chronological order &,;, i € 
{1,2,...,switen} with the total number of switches denoted by Nswitch. The time index 
given by the first sampling point of the measurement data represents the first occurrence of 
the initial state and is thus defined to be k, ı = 1. 
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The time index right before a switch, i.e. the last time index the current state is active, is 
called end of the current segment kj; i = kr,i+1 — 1. The end of the last segment is given by 
the last time index of the measurement data T which leads to k = T. A schematic 
sketch for Nswitch = 3 is given in Fig. 6.2. 


y. 
T ;Nswitch 


Segment 1 Segment 2 Segment 3 


ka kri ko kr 2 ka ka 


Figure 6.2: Schematic view of switches 


6.1 Verification of Hybrid Systems with Mapped State 
Signal 


The introduced definitions can be used to set up the specification of a hybrid system. 


Definition 6.8 (Specification of an Interval Type Hybrid System) 
An interval type specification of a hybrid system according to Def. 6.6 is given by 


St, :— (25,57) (6.13) 


with the nominal state machine Z* according to Def. 6.4 and the set of nominal dynamic 
systems S*, according to Def. 5.1. 


A special case of Def. 6.8 is given if the parameters of the dynamic subsystems are point real 
values. 


Definition 6.9 (Specification of a Point Real Hybrid System) 
A point real specification of a hybrid system according to Def: 6.6 is given by 


S5 4:7 (25,83) (6.14) 


with the nominal state machine Z* according to Def. 6.4 and the set of nominal dynamic 
systems 57, according to Def. 4.2. 


The dynamic part of a system according to Def. 6.9 is also called linear time variant system 
with piecewise constant parameters. 
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The verification method developed in this chapter is again based on interval enclosures of 
the input and output measurement data 


\Uneass Yineas| = [cee m , ees m (6.15) 


extended by a measured state signal Qmeas according to Assumption 6.1. The state signal 
Qmeas leads to a set of measured states Qmeas. It is assumed that the elements of the set of 
measured states can be mapped on the set of specified states. In this case, the set is called 
mapped set of states according to Def. 6.10. 


Definition 6.10 (Mapped Set of States) 
The set of measured states Qmeas, consisting of the unique values of the state signal 
Qineas = Cae ee is called mapped with the specification S7, ,, if all measured states 


ados € Omeas with j € {1,2,...,|Qmeas|} can be mapped to an equivalent nominal 


state q9))* € Q* given in the specification i.e. 


as = ghar, J € 11; 2, ne |Qmeas|}- (6.16) 


As S% q is a special case of SH Def. 6.10 is also valid for SH. In case Qmeas is a mapped 
set of states, the state signal Qmeas is called mapped state signal. If the mapped state q')* 
is known, the respective mapped nominal subsystem s'(/)* is also known. The hybrid verifi- 
cation problem can now be formulated as follows: 


Problem 6.1 (Mapped Set of States Based Point Real Hybrid Consistency) 
Is the nominal hybrid system, specified by a point real hybrid specification 


SHa = 12,54} (6.17) 


consistent with the input-output behavior given by the interval type enclosures of T' mea- 
surement values 


HA T 
[U neas Yıneas] = s , (Ymeas,k) p1 (6.18) 


and the measured mapped state signal Qmeas, i.e. can the measurement data be explained 
by the nominal system? 


The problem can be solved by tackling the two system parts individually. This is possible due 
to an implicit connection given by the matching of the discrete states which is done based 
on the dynamic parameters. First the dynamic subsystem is considered and verified in the 
next section. Afterwards the verification of the discrete event part is introduced. Finally the 
results are combined to verify the overall hybrid system. 
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6.1.1 Verification of the Dynamic Subsystems 


To verify the individual dynamic subsystems it is necessary to split the interval type mea- 
surement data [U meas, Ymeas| based on the information given in the mapped state signal. 
The resulting segments 


; T Nswitch k; j 
( bm Sal ) = ( en 2a 


ks : Nswitch 
T RE (Ymeas,k) kk, 3l ) (6.19) 
I= j 


j j=l 
can be verified individually against the respective specification in the set of subsystems 57. 
It is possible to verify the individual dynamic behavior? of each state present in the measure- 


ment data qus € Qmeas as defined in Prop. 6.1. 


Proposition 6.1 (Dynamic Consistency of a Segment) 

The interval type enclosure of the measurement data is split into j € {1,2,... ‚Ngwiten} 
parts, given by the segments lom YO). ; Each segment repre- 
sents a specific state qlas, a regressor matrix AQ; and a  measure- 


ment vector Bias. The segment j is dynamic consistent with the respec- 
tive mapped subsystem s“O))* € à df the specified system parameters 


j ; . . * . . * . . * . * * . y * . . T 
eC» = lai (i(3)). a» (i(j)), EE: ET) Cy ((3)). EXER EET Co a EU) 
are part of the united solution set y = (A ta Bas), Le. if 


QG»* € Y (6.20) 


Proof: 
The measured state signal Q,ncas includes the true active states. It is thus guaranteed that 


only measurement data generated by subsystem sS. based on the active state glas is 


included in uT YG s and that this data is not corrupted by measurement generated 


(i) 


by other subsystems. The current state gr... is a mapped state according to Def. 6.10 and 
thus the connection between measurement and specification is also correct and the respective 
nominal subsystem 5 ))* is known. 

Using this information, the setting can be reduced to the time invariant consistency problem 
given in Problem 4.1 for segment j and subsystem s /))* and time invariant consistency can 
be checked according to Prop. 4.1 which proves Prop. 6.1. 


The considerations are now extended to the complete set of measurement data, consisting of 
several segments. 


* As a direct specification 7 is used in the definition, “dynamic consistency” means “time invariant full consis- 


tency” throughout this chapter. 
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Proposition 6.2 (Dynamic Consistency of the Measurement Data) 
The mapped state signal Qmeas and the segmented measurement data 
; 5 Nswitch ; 
( [Ufa X ) of a VO, leading to the state las, the regressor matri- 
j=1 

ces AU) s and the measurement vectors BY... with j = {1,2,...,Nswiten}, are 
dynamic consistent with a set of direct specifications Sž, if there is dynamic consistency of 
each segment given in the measurement data with its respective mapped subsystem s(U)* 
Le: 


EM e 3^0). vi € (1,2,... nswach]- (6.21) 


Proof: 
Prop. 6.2 results straight forward by applying Prop. 6.1 to all j = (1,2,..., nswitcn] seg- 


Nswitch 


ments given in the segmented measurement data ( Be ) 


je 


It is possible to reformulate this proposition to get the inverse relation similar to Prop. 5.2. 


Proposition 6.3 (Inverse of Dynamic Consistency of the Measurement Data) 
The mapped state signal Qmeas and the segmented measurement data 


(i) G) REDEN ; G) = 
Urmeas, Y meas| ) . of a VO, leading to the state qmeas, the regressor matri 


ces AG). and the measurement vectors B®... with j = {1,2,...,switen} are called 
dynamic inconsistent with a set of direct specifications Sž, if there is at least one segment 
j € {1,2,...,switen} given in the measurement data that does not show dynamic 
consistency with its respective mapped subsystems s * according to Prop. 6.1, i.e. 


LLI 


$12, inguina] | OOO g SS (6.22) 


Proof: 
For dynamic consistency of the measurement data according to Prop. 6.2 it is necessary that 
all segments of the measurement data are dynamic consistent according to Prop. 6.1. If there 
is a segment that does not show dynamic consistency, this segment can not be explained 
by the specification. Thus there is unspecified behavior and it is not possible to explain the 
whole measurement data by the specification. 


Note that these propositions are based on the segments given in the measurement data. How- 
ever it is possible that there is nominal behavior that is not present in the measurement data. 
Though this will not change the dynamic consistency result for the measurement, it will 
influence the discrete event verification result introduced in the next section. 
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Also note that a mapped state signal is used in Prop. 6.2 and Prop. 6.3. This means that all 
segments given in the measurement data can be mapped to the specification, as stated in 
Def. 6.10. 


6.1.2 Verification of the Discrete Event System 


The next step in the verification of the hybrid system is given by regarding the discrete event 
part. Therefore it is necessary to determine the state machine Zmeas that generated the 
measurement data. The active states given in the mapped state signal Qmeas represent a trace 
of the unknown generating state machine Zmeas. The system dynamics represent the state 
of the system and not an emission of a state or an event. Therefore it is in this case possible 
to reconstruct Zmeas based on this trace. The reconstructed generating state machine Zneas 
can then be compared with the specified state machine Z* by comparing the corresponding 
states and transitions. 


Proposition 6.4 (Full State Consistency) 


The mapped state signal Qmeas given for the discrete sampling points k = {1,2,...,T}, 
leading to the mapped set of states Omeas of the discrete event part Zmeas of a VO is called 
full state consistent with the nominal state machine Z*, if 


Q* = O ineas: (6.23) 


Proof: 
All nominal states Q* are given in the specification of the state machine Z*. The states 
implemented in Zmeas are given in the mapped set of states Omeas. According to Def. 6.10, 
this means that both sets are defined on the same elements. Full state consistency means that 
exactly the specified states are given in the measurement data, i.e. that both sets contain the 
same elements. This comparison is given in (6.23) and proves full state consistency. 


Due to a measurement scenario that does not cover all states it is possible that not all dynam- 
ics are present in the measurement data. This case leads to partial state consistency. 


Proposition 6.5 (Partial State Consistency) 


The mapped state signal Qmeas given for the discrete sampling points k = {1,2,...,T}, 
leading to the mapped set of states Omeas, of the discrete event part Zmeas of a VO is called 
partial state consistent with the nominal state machine Z*, if 


Q" D Omeas- (6.24) 
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Proof: 
Based on Prop. 6.4 and Def. 6.10, Qmeas and Q* are defined on the same elements, i.e. there 


is a mapping as = qË for j € {1,2,...,|Qmeas|}. If additionally (6.24) holds, 


adus € Qmeas > aas = gi e Q“ (6.25) 


holds as well. This means that all implemented states qus € Omeas are also part of the 
specification Q*. Thus there are only specified states in the measurement data. However, 
there are states in the specification that are not part of the implementation and prevent full 
state consistency. The discrete event part Zmeas of a VO is hence called partial state consis- 
tent. 


Partial consistency of the discrete event part Zmeas of a VO is a hint to improve the measure- 
ment scenario or to collect more measurement data. 

If there is neither full nor partial consistency, the system is state inconsistent according to 
Prop. 6.6. 


Proposition 6.6 (State Inconsistency) 


The state signal Q meas given for the discrete sampling points k = {1,2,...,T}, leading to 
the set of states Qmeas, of the discrete event part Zmeas of a VO is called state inconsistent 
with the nominal state machine Z*, if 


Q* Z Omeas: (6.26) 


Proof: 
If (6.26) holds, there are implemented states glas € Qmeas that are not part of the speci- 
fication O*. Thus there are unspecified states in the measurement data. The discrete event 
part Zmeas of a VO is hence called state inconsistent. 


Prop. 6.6 can be connected with the inverse of dynamic consistency of the measurement data 
as given in Prop. 6.3. If there is an additional state "M Q*, the respective measurement 
data lum IM cannot be explained by any s* within the specification. In this case 


there is both, dynamic inconsistency according to Prop. 6.3 and state inconsistency according 
to Prop. 6.6. 


The second part of the discrete event system to be verified is the transition function. Addition- 
ally to the mapped state signal’ a the set of measured events Emeas needs to be 


G) € Emeas if Wmeas,k = DR 


obtained. According to Def. 6.2 there is an active event e,, ..; d 


? The last measurement value for k = T can not be evaluated as there is no following state k = T + 1. 
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Full transition consistency is then defined as follows: 


Proposition 6.7 (Full Transition Consistency) 


The mapped state signal Ces dias and the set of events Emeas of the discrete event part 


Zmeas Of a VO are called full transition consistent with the nominal state machine Z*, if 


a 


(Vameas,s € (amens k) a ) ^ (3 N sk = — : du hasi = qmeas,k+1)- 
(6.27) 


Proof: 
All nominal transitions are given in the transition function f* according to Def. 6.3. The tran- 


sition function f*(qmeas,k, e a „) is evaluated for the current measurement state ¢meas,k 


and the current events e? 


eneas k- É the transition function yields the following measurement 


state qmeas,k+1 for at least one event eO x the right hand side of (6.27) holds. Thus the 
observed transition at time k is part of the nominal transition function. 
If the right hand side of (6.27) holds for the measurement sequence k = (1,2,..., T — 1}, ie. 


(Ydmeas,k € (dmeas, PELi), all observed state transitions are defined in the nominal transi- 


tion function f*. Thus the measurement is full transition consistent. 


Prop. 6.7 also implies that the current values of the enabler signal Wmeas,, are within the 
nominal limits [* at each switch k = kr; — 1 with i = {2,3,...,switcn}.!° Other than 
state consistency, full transition consistency can be achieved although a nominal transition 
is not triggered by the measurement data. 


The notion of partial transition consistency includes specified transitions that are triggered 
at unexpected times. 


Proposition 6.8 (Partial Transition Consistency) 


The mapped state signal DNE: da and the set of measured events Emeas of the discrete 


event part Zmeas of a VO are called partial transition consistent with the nominal state 
machine Z* , if 


| 


Qmeas,k € er N (36 = 2) : OF eau, é) = Üosus edt.) ^ (é 4 ee) 
(6.28) 


10 It is not necessary to check the first switch, as it represents the begin of the experiment. 
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Proof: 
The transition function f* is evaluated as explained in Prop. 6.7, but using the set of nominal 
events £* instead of the set of measured events Emeas. Thus f*(qmeas.k, €) can yield the 


correct following measurement state qmeas,k+1 for any specified event € € €*, regardless of 
(i) 

meas,k* 

Condition (6.28) is fulfilled if there is at least one transition in the measurement sequence 
(i) 


meas,k* 


the current measured events e 


k = (1,2,..., T — 1} that was triggered by an unexpected event č 4 e 


Unspecified transitions and transitions connecting unspecified states lead to transition in- 
consistency according to Prop. 6.9. This represents the situation, where it is not possible to 
explain the observed transition by the transition function. 


Proposition 6.9 (Transition Inconsistency) 


The mapped state signal Co ee and the set of events Emeas of the discrete event 


part Zmeas of a VO are called transition inconsistent with the nominal state machine Z*, 


if 


(Jimenzk € a) N (ve € £*) : (f* (dmeas ki é) # RE (6.29) 


Proof: 
Condition (6.29) is fulfilled if there is at least one unspecified transition at a state qmeas,k € 


RN ey that is not defined for any nominal event € € £*. Therefore the measured 
transition is unspecified, which is inconsistent in the sense of this thesis. 


The results of state and transition consistency are combined such that the weakest result 
of the individual conditions determines the result of the whole discrete event system. The 
respective propositions are given in the following. 


Proposition 6.10 (Full Discrete Consistency) 


The discrete event part Zmeas of a VO is called full discrete consistent with the nominal 
state machine Z*, if there is full state consistency according to Prop. 6.4 and full transition 
consistency according to Prop. 6.7. 
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Proof: 
The main components of the discrete event part Zmeas according to Def. 6.4 are the set of 
states Q*, the finite set of events €* and the transition function f*. The consistency of 
this components is verified using the propositions about state and transition consistency. 
Only if all components show full consistency according to the respective propositions, the 
overall system also shows full consistency. Thus full discrete consistency is only given if 
there is full state consistency according to Prop. 6.4 and full transition consistency according 
to Prop. 6.7. 


Proposition 6.11 (Discrete Inconsistency) 


The discrete event part Zmeas of a VO is called discrete inconsistent with the nominal state 
machine Z*, if there is state inconsistency according to Prop. 6.6 or transition inconsistency 
according to Prop. 6.9. 


Proof: 
Discrete inconsistency is also an integral property of the main components of the discrete 
event part Zmeas as shown in Prop. 6.10. Therefore the overall discrete system is inconsis- 
tent as soon as there is at least one inconsistent component, i.e. if there is state inconsistency 
according to Prop. 6.6 or transition inconsistency according to Prop. 6.9. 


Proposition 6.12 (Partial Discrete Consistency) 


The discrete event part Zmeas of a VO is called partial discrete consistent with the nominal 
state machine Z*, if there is neither full discrete consistency according to Prop. 6.10 nor 
discrete inconsistency according to Prop. 6.11. 


Proof: 
If there is neither full discrete consistency according to Prop. 6.10 nor discrete inconsistency 
according to Prop. 6.11, there is at least one main component that shows partial consistency 
according to Prop. 6.5 or Prop. 6.8. Even though the other main component might show full 
consistency according to Prop. 6.4 or Prop. 6.7, the integral property can not be better than 
the properties of the included main components. 


The implementation of this propositions is straight forward, as all necessary sets and values 
are available in the used setting. More realistic scenarios assuming less a priori knowledge 
are introduced in Section 6.2 and Section 6.3. 
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6.1.3 Combination of the Dynamic and the Discrete Verification 
Results 


The results of the discrete event system part can be joined with the results of the dynamic sys- 
tem part to achieve the overall assessment of the hybrid system. The solution of Problem 6.1 
is thus given by Prop. 6.13 as follows: 


Proposition 6.13 (Mapped Set of States Based Point Real Hybrid Consistency) 


The mapped state signal Qmeas, the interval type enclosures of the measurement data 
[(Umeas, Ymeas| and the set of events Emeas of the verification object Hmeas are called 
consistent with a direct hybrid specification S7, 4, if 


* the specified state machine Z* is full discrete consistent and 


e the specified subsystems s\* € 5; are full dynamic consistent. 


Proof: 
Based on the mapped state signal according to Def. 6.10 the currently active subsystem is 
known at each time step. The resulting trace of the state machine is checked for consistency 
with the discrete event specification Z* by regarding the states and transitions. Consistency 
of the states can be checked by Prop. 6.4 and consistency of the transitions by Prop. 6.7. If 
both parts are consistent, the measurement is consistent with the specified state machine. 

The dynamic of the included subsystems can be checked according to Prop. 6.2, again assum- 
ing a mapped set of measured states Qmeas. If both system parts are consistent, the overall 
hybrid system is consistent. 


Consistency of the hybrid system results from the combination of results for the dynamic and 
discrete subsystems according to Prop. 6.13. Inconsistency of the hybrid system is given as 
soon as one subsystem is inconsistent. Else the hybrid system is partial consistent, meaning 
that there is no inconsistent subsystem but at least on subsystem is partial consistent. 
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An overview of the different propositions and possible results is given in Tab. 6.1. 


Table 6.1: Consistency criteria 


Property 

Part of the System full consistency | partial consistency | inconsistency 
Dynamic Subsystems Prop. 6.2 not applicable Prop. 6.3 
States Prop. 6.4 Prop. 6.5 Prop. 6.6 
Transitions Prop. 6.7 Prop. 6.8 Prop. 6.9 
State Machine Prop. 6.10 Prop. 6.12 Prop. 6.11 

Prop. 6.3 
Hybrid System Prop. 6.13 else or 

Prop. 6.11 


By applying Prop. 6.2, consistency of the dynamic part of the hybrid system can be shown in 
a straight forward way. The results are calculated based on the propositions of the previous 
chapters and thus show the same guaranteed properties as defined there. A direct specifi- 
cation Sž was used throughout the chapter for notational simplicity. It is straight forward 
to extend all propositions to hold also for an interval type specification 57. This is due to 
the fact that the dynamic verification is based on the united solution set given by the mea- 
surement data. Therefore all properties stay the same except that the verification of dynamic 
consistency is done based on Prop. 5.1 instead of Prop. 4.1. 

Subsystems that are verified using the Kaucher based approach are guaranteed to be correct 
and it is not possible that there are any hidden faults present in the system. 

The state machine has to use point real numbers in both cases. Therefore the definitions and 
propositions for the discrete event part are unchanged. 

An example of the hybrid verification procedure for a mapped set of measured states is given 
in Example 6.1. 
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Example 6.1: 


Assume a direct hybrid specification S% , = [Z*, S7]. The nominal state machine Z* is 
depicted in Fig. 6.3. 
3) 


4) 


el el 


) 


e? 

20D O 
NT 

Figure 6.3: Specified state machine Z* 


The state machine consists of two states Q* = (q?*, ar} = {1,2}, the events are based 
on the enabler signal that is defined to be the output signal w = y: 


e(0* : wp ei = n, 2] (6.30) 
e?)* : wg € 12) = [31, 33] (6.31) 
e(9* : wp € I8) = [-oo, oo] (6.32) 
e(9* : wp € UY = [-oo, oo]. (6.33) 


The events e(®)* and e(9* are enabled for all values of the enabler signal, leading to the 
permanent possibility to stay in the current state, based on the transition function: 


f (a, ee) = 2 (6.34) 
F (a, e) = gi) (6.35) 
f (a, e) = gO (6.36) 
f (2, eC) = q. (6.37) 


The set of dynamic systems S is given by first order systems, ie. Nali) = ne(i) = 1, 
Vie {1,2} leading to 


s(* > Yk = ay (i) yR—1 + ci (i)Up_-1.- (6.38) 


The nominal parameters of the dynamic subsystems are given in Tab. 6.2. 


Table 6.2: Nominal parameters of the subsystems s‘!)* and s(2)* 
Subsystem | aj | cT 


qU* 0.1| 1 
qO* 20 | 1 
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There are no optional initial input and output values given. Thus the input and output values 
are kept across the switches. 

The implementation is assumed to be done by one or more human developers. Therefore there 
might be inconsistencies in the resulting VO. Note that the implemented system is assumed 
to consist of real hard- and software and to include a given plant that cannot be changed. 
Therefore the corresponding state machine Zmcas and its dynamical subsystems Smeas are 
not directly known. Nevertheless it is possible to excite the system and measure its output and 
state signals. The random excitation signal used in this example is given by 


uj = 1 + 0.2%, (6.39) 
where 1j, is drawn from a standard normal distribution. The resulting measurement data 


is given in Fig. 6.4. Thereby the output was enclosed by intervals using an additive fault of 
dy, = 0.5. The switching times are based on the information given in the mapped state signal 


Q meas = (1, 1, 1, L 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1] (6.40) 
leading to the switches k- = [1,6, 10]. The relevant values of the enabler signal are given 
at the time steps right before a switch wy, ,-1 withi € (2,3) ie ws = [0.6, 1.1] and 


wo = [31.9, 32.9). 


2 


| | | 
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 


gj Erts mice toe ----- |—— Measurement |- | 
—— Enclosure 


1 2 3 4 5 6 7 8 9 10 1 12 13 M 15 
Time step k with At — 1s 


Figure 6.4: Measured trajectory and subsystem switches 
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Dynamic Consistency 


The subsystems are verified using the introduced Kaucher based method. The resulting fea- 
sibility signals are depicted in Fig. 6.5. The first verification result can be calculated for 
kmin = max(na,ne) + 1. Due to the autoregressive system of order n, = Ne = Lin 
this example it is thus not possible to calculate a verification result for the very first element 
of each segment. It can be seen, that all three segments j = {1, 2,3} can be explained by the 
respective mapped nominal states i = {1,2} using the Kaucher based method according to 
Prop. 6.2. Therefore the dynamic subsystems Seas are full consistent with the specification 
Si. 


T T T T T T T T T 
Verification Result Segment 1 and s 


feasible -- --m---m---m---g---- N 
BEN 8 
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 


T T T T T T T T T T 
Verification Result Segment 2 and s?* 
feasible 


MM N 
infeasible NN 
1 2 3 4 5 6 v4 8 9 10 11 12 13 14 15 


T T T T T T T T T 
Verification Result Segment 3 and s 


nn mn. m---m--------# 


feasible 
infeasible Ammann | 


|? i i 
1 2 3 4 5 6 7 S 9 10 11 12 13 14 15 


Figure 6.5: Verification result for each segment (using a mapped state signal) 
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Discrete Consistency 


There is a mapped set of states given in this example. Therefore the measured set of states is 
given by the unique values in the mapped state signal Qmeas 


One 12) = fa ou que = Q* (6.41) 


which leads to full state consistency according to Prop. 6.4. 
The set of measured events can be extracted from the measurement data. First, the events 
leading to the switches in k € {6, 10} are determined: 


e(D s: ws = (0.6, 1.6] n [1 2] 40: e(?* (6.42) 

ficus, : Ws = [0.6, 1.6] N [-00, oc] # 0 : c9? (6.43) 

eG) 5: Ws = (0.6, 1.6] n [-oo, oo] Z 0 : e(?* (6.44) 

e(D so: wo = [31.9, 32.9] n [31, 33] 40: e®* (6.45) 

as: Wo = [31.9, 32.9] n [7co, oo] #0: e®* (6.46) 

as: Wo = [31.9, 32.9] n [7co, oo] ZU: e(?* . (6.47) 

This means that etus. eG», ex} are activated atk = 5 and [e e(3)* ,e(4)* } are 
activated at k = 9. The events can now be applied to the nominal iransition ee 

* (roae an) =f* (a ! *, e) + q?* = dmeas,6 (6.48) 


2 * ! * 
(Qmeas,5; en) NN = f (4 1 i eO =q 0) # dmeas,6 (6.49) 


) 

qt eo) = 0 (6.50) 
) = Gmeas,i0 (6.51) 
) 


( 
Qmeas,9: €meas,9 


q 
= (6.52) 


(3) ) 
(Qmeas,5; Cmeas,5 


! 
q d r elt») = q # Gmeas,10- (6.53) 


It can be seen that (6.48) and (6.51) hold. This leads to the verification of the transition 
function in k € {5,9}, that generates the switches k, € {6,10}. The events e®)* and e4)* 
are enabled for allk € {1,2,...,15} and provide the possibility to stay in the same state 
for several time steps. To improve readability, only the verification of one exemplary step is 
shown. The activated events are 


ee 3 i Ws = [15.3, 16.3] n [—co, oo] #9: a (6.54) 
eC) gg: ws = (15.3, 16.3] n [-00, oo] £0: e®*. (6.55) 


These are applied to the nominal transition function 
J* (qmeasi8iCmonne) =f" (a, =0 (6.56) 


* 2 * * * ! * 
f C =f C , e ) = qO = dmeas,9- (6.57) 
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Thus the behavior is valid at time k = 8. Similar results are obtained for all other time steps. 
The nominal transition function holds Vdmeas,k € Qmeas, which leads to full transition 
consistency according to Prop. 6.7. Hence full state consistency (Prop. 6.4) and full transi- 
tion consistency (Prop. 6.7) hold in the given example. This leads to full discrete consistency 
between Zmeas and Z* according to Prop. 6.10. 


Hybrid Consistency 


The previous partial results are now combined with respect to hybrid consistency as given in 
Prop. 6.13. It was shown that Zmeas is full consistent with Z* and that Smeas is full consistent 
with S3. Therefore the verification object Hmeas and the specification S7, q are consistent. 
This means that the superimposed state machine as well as the linear dyna mic subsystems of 
the VO that produced the measurement in Fig. 6.4 are full consistent with the specification of 
Si q given in this example. Therefore the VO is verified with respect to the nominal system. 
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6.2 Verification of Hybrid Systems With Given Switching 
Times 


In general it is not possible to measure the internal signals and states of a VO. Therefore the 
setting is changed and the assumption of an available mapped state signal Q meas is dropped. 
Nevertheless, it is still assumed that the correct times of the switches k, are available, even 
though the respective active states are unknown. The resulting consistency problem is given 
in the following: 


Problem 6.2 (Point Real Hybrid Consistency with Given Switches) 
Is the nominal hybrid system, specified by a direct hybrid specification 


Sha — 125.81) (6.58) 


consistent with the input-output behavior given by the interval type enclosures of T' mea- 
surement values 


T T 
[E eaa; Y ricas] = [mut , (Ymeas,k)k—1 (6.59) 
and the set of switches 
{keg heey (6.60) 


ie. can the measurement data be explained by the nominal system? 


To solve the problem, it is necessary to determine a mapped state signal Q. The measurement 
data can be segmented based on the given switches. The result can be interpreted as an 
unmapped state signal as it is not known which state is active after each switch. Therefore 
is is necessary to determine the correct nominal subsystem for each segment. There are two 
important assumptions that need to hold to ensure an unambiguous mapping from nominal 
subsystems to measurement segments. 


Assumption 6.2 (Prager-Oettli-Distinguishability of the Set of Subsystems) 


Each two nominal subsystems s*, s(* € Sj with i # j are called Prager-Oettli- 

Distinguishable (PO-Distinguishable) with respect to specific segmented measurement data 

ze Ye) nn if there is no segment j € {1,2,...,Ngwiten} that fulfills 
j=l 

Prop. 4.1 for both subsystems s(?* and s9)*, 


This means that two distinct nominal subsystems given in the specification are sufficiently 
different with respect to the measurement data, noise assumptions and the interval enclosure. 
The second assumption transfers this property to the segments of the measurement data. 
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Assumption 6.3 (Mappability of the Measurement) 


The segmented measurement data ( fusa | ) ns is called mappable to the set 
j=l 

of dynamic subsystems S*, if there is no specified subsystem s®* € S; that fulfills Prop. 4.1 

for any segment Be aT that was generated by another dynamic system s) with 


PORO 


If Assumption 6.2 and 6.3 hold, it is possible to determine a mapped state signal. Therefore 
Problem 6.2 can be solved by the following proposition: 


Proposition 6.14 (Point Real Hybrid Consistency with Given Switches) 

The direct hybrid specification S7, 4, the interval type enclosures of the measurement data 
[U neas; Ymeas| and the set of switches {kes} of a verification object Hmeas can 
be used to set up a mapped state signal Q if the specified set of subsystems S7 is PO- 
distinguishable with respect to the measurement data according to Assumption 6.2 and the 
measurement data is mappable to S} according to Assumption 6.3. 

The availability of a mapped state signal Q transforms the problem to a mapped set of states 
based point real hybrid consistency problem. This problem can be solved using Prop. 6.13. 


Proof: 
An unmapped state signal can directly be constructed from the given correct switches 
{krj par” according to (6.19). As Assumption 6.2 and 6.3 hold, individual nominal sub- 
systems within the specification Sž can be distinguished from each other. 

Also, the united solution set Y defined by the measurement data of each segment 
prm YO). , generated by a dynamic system sU). cannot be explained by any other 
nominal system s“* € Sž, with sU) 7 s(?*, The generating subsystem can thus be de- 
termined unambiguously for each segment. The mapped nominal subsystems can be used 
to determine the active states q, with k € {kr j, krj +1,...,k, j} for all segments j € 
{1,2,...,switen}. The resulting signal Q is a mapped state signal according to Def. 6.10. 
The given measurement data [Umeas, Ymeas|, together with the given set of nominal subsys- 
tems S* and the extracted mapped state signal Q represents the setting of Problem 6.1 that 
can be solved by Prop. 6.13. 
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Segments generated by an unspecified subsystem 5 ¢ S7 cannot be mapped to a nominal 
subsystem s(?* if § is PO-Distinguishable from all nominal subsystems s()* € Sj according 
to Assumption 6.2. If this is not the case, the measurement data generated by 5 can be ex- 
plained by at least one nominal subsystem s(?* € 5; and thus it is impossible to recognize 
the subsystem 3. Note that the definition is based on all measurement data of a segment, 
ie. it is possible that partial measurement data of a segment can be included in more then 
one nominal subsystem. 


The method leads to a mapping algorithm that compares the dynamic of each nominal subsys- 
tem, given by the n, states, with every segment j € (1,2,..., Nswitch } of the measurement 
data. This comparison is done based on the united solution set, as given in Prop. 4.1. The 
flowchart of the algorithm is given in Fig. 6.6. 

An exemplary application of point real hybrid consistency with known switches is given in 
Example 6.2. 


No matching state 
for segment j 


Assign state i 
to segment 7 


Figure 6.6: Flowchart of the mapping algorithm 
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Example 6.2: 


Assume the same setting as introduced in Example 6.1, without the assumption of a given 
measured mapped state signal. Instead, a set of switches is given by 


k- € {1,6,10}. (6.61) 


Based on the switches and the length T = 15 of the measurement data, the endpoints of the 
segments can be calculated: 


kr € {5,9,15}. (6.62) 


The nominal parameters given in Tab. 6.2 and the algorithm of Fig. 6.6 are used to determine 
the mapping between segments and states. 

To show the correct classification of all subsystems, the algorithm is altered such that it com- 
pares all possible nominal subsystems with each segment instead of proceeding to the next 
segment as soon as a consistent subsystem is found. 

The results are depicted in Fig. 6.7. Each subfigure shows the evaluation of both nominal 
subsystems s“)* and s'2)* for one of the segments j € {1,2,3}. The shaded areas mark 
measurement data that does not belong to the segments and thus is not taken into account for 
the respective verification. It can be seen, that only one generating nominal subsystem can 
be verified for each of the three segments. This leads to an unambiguous mapping 


Q-—]11,1,1,1,9,2,3.2. 11,111; 1] (6.63) 


which is the same as given in (6.40). 

The matching algorithm is based on Prop. 6.1. Therefore dynamic consistency is given, if 
it is possible to generate a mapped state signal Q. The discrete verification can be done as 
described in the previous section. The combination of both yields the hybrid verification result. 
Both results are equivalent with the calculations and results demonstrated in Example 6.1 and 
are thus not repeated here. 
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Figure 6.7: Verification Result for each segment (without mapped set of states) 


6.3 Verification of Hybrid Systems With Unknown 
Switching Times 


In the last step even the knowledge about the switching times is dropped. Only the specifi- 
cation and the input and output measurement data are available. Thus the switching times 
as well as the respective modes need to be reconstructed from the measurement data only. 
Therefore an additional segmentation step is necessary in the procedure. This step aims at 
finding the unknown switches k+, at which the active generating subsystem changes. The 
corresponding verification problem is given in Problem 6.3. 


6.3 Verification of Hybrid Systems With Unknown Switching Times 
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Problem 6.3 (Point Real Hybrid Consistency) 
Is the nominal hybrid system, specified by a direct hybrid specification 


Sia = 12,54} (6.64) 


consistent with the input-output behavior given by the interval type enclosures of T' mea- 
surement values 


T F 
|U meas Yıneas] = ra , (Umeasiki r (6.65) 


without any knowledge about the state or the switching times? 


The problem can be solved, if it is possible to determine a mapped state signal from the input- 
output data only. If this signal is constructed, Problem 6.3 can again be reduced to Problem 6.1 
which is solved by Prop. 6.13. However, in this setting there is no information about the active 
states or the switches. Therefore an additional identification and segmentation procedure is 
necessary. The procedure introduced in this thesis is based on previous work of the author, 
published in [Die13a], [Die13b] and used in [Die17]. However, while the original work aims 
on identifying an a priori unknown library of subsystems from the measurement data, the 
concept is adapted in this thesis to consider a given set of nominal subsystems. Additionally, 


it is extended to the interval arithmetic context of guaranteed verification. 
First, switches that are detected by the segmentation and identification method are defined. 


Definition 6.11 (Detected Switch) 


The time instant k that does not show dynamic consistency according to Prop. 6.1 of any 
nominal subsystems s(0* € S3 with the current regressor matrix Ameas,k and the current 


measurement vector Bmeas,k is called detected switch k, ;+1 = k. 


To ensure correct segmentation, the available measurement data needs to be Prager-Oettli- 


Segmentable: 


Assumption 6.4 (Prager-Oettli-Segmentability) 
The measurement data Oe ads Sa ae and luccm of two consecutive seg- 


ments j and j + 1 that are generated by distinct subsystems sU)* + sUÜ+V* are called 
Prager-Oettli-Segmentable with respect to a given set of subsystems $7, if there is no sub- 
system s()* € 5; that fulfills Prop. 6.1 for all measurement data points in the combined 
segment |kr j, kr j+1]- 


Assumption 6.4 implies Prager-Oettli-Distinguishability of the set of subsystems S% accord- 


ing to Assumption 6.2. Based on the segmentability assumption it is possible to extract 
mapped state signal to solve Problem 6.3. 


a 
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Proposition 6.15 (Point Real Hybrid Consistency) 


The direct hybrid specification S7, q and the interval type enclosures of the measurement 
data [Umeas,; Ymeas| of a verification object Hmeas can be used to set up a mapped state 
signal Q, if the measurement data is PO-segmentable with respect to S7, q. 

The availability of a mapped state signal Q transforms the problem to a mapped set of states 
based point real hybrid consistency problem. This problem can be solved using Prop. 6.13. 


Proof: 

If Prager-Oettli-Distinguishability (Assumption 6.2) and mappability of the measurement (As- 
sumption 6.3) are fulfilled, there is only one possible consistent subsystem s)* € S 7 at the 
end of a segment kept = kr jti — 1. This subsystem needs to represent the active state of 
the whole segment es kr : 

On the other hand, it is impossible that this subsystem is consistent with the measurement 
data of the entire next segment generated by a different subsystem s@+1), A nominal sub- 
system s(?* € Sï that fulfills Prop. 6.1 for all measurement points in the entire combined 
segment [k- j, kr’ 341] will also fulfill Prop. 6.1 for any part of the combined segment. This 
holds especially for the parts [k+ j, krj] and [kr ;+1, kr’ j+1], ie. the same nominal subsys- 


tem is consistent with two distinct segments [Ua er CUM and Br ; yor) . This 


is only possible if s®* = s@) = s@+1) which is not allowed in the case of mappability of 
the measurement according to Assumption 6.3. 

Thus it is possible to determine a switch kr > kz j for each j € {2,3,..., Nswiten}, Le. the 
true number of segments is determined even if the detected switches do not exactly match the 
real ones. Based in this result it is possible to determine the respective active subsystem for 


all segments, whereas each estimated segment |i. js ks. ] at least partly overlaps with the 


true segment [k+ ;, krj]. This leads to a mapped state signal with correct mapping for all 
segments, even if the detected segment boundaries might slightly differ from the true ones.!! 
Thus Problem 6.3 is transformed to Problem 6.1 which completes the proof. 


If Prager-Oettli-Segmentability according to Assumption 6.4 does not hold, it is not possi- 
ble to determine the switch. The flowchart of the algorithm implementing the introduced 
identification and segmentation method is given in Fig. 6.8. 


11 A detailed proof of this property is given in Section 6.3.1. 
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Figure 6.8: Flowchart of the identification and segmentation algorithm 


The inner loop (’i-loop’) depicted in the flow chart compares all n; = |S}| nominal subsys- 
tems s(?* € S*, with the current segment. If there is a consistent nominal subsystem, it is 
added to the set of consistent subsystems Scon. The next loop (’k-loop’) is running as long as 
the set of consistent subsystems Scon,k is non-empty. If this is the case, the current segment 
can be explained by at least one nominal subsystem. Therefore the segment is extended by 
one time step. 

This new segment is again verified by the i-loop. If the set of consistent subsystems is empty, 
ie. Scon = 0), none of the nominal subsystems is able to explain the current segment. How- 
ever, the measurement was verified in the previous step. Therefore a detected switch’ is 
recognized at krja = k and the active subsystems for the segment j are given by Scon,k-ı- 
Note that due to Assumption 6.3 only one subsystem is allowed to explain an entire segment 
L1 

However, multiple consistent subsystems |Scon,k| > 1 are possible for partial segments 
ke {krj +1, krj +2,..., krj = 1}. 

Detecting a switch and determining the state of the finished segment ends the k-loop of the 
algorithm in Fig. 6.8. The measurement values belonging to the just finished segment are 
removed from the considered measurement data. 


of the measurement data. This leads to 5 A 
con, ky’ j 


12 The first switch of a system is defined to be k,,1 = 1, according to Def. 6.7ff. 
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All counters and intermediate values are reinitialized and the next iteration of the outer loop 
(j-loop’) starts for the following segment. 


In order to achieve the hybrid verification result, the discrete and dynamic results are com- 
bined according to Section 6.1.3. An application of the procedure is given in Example 6.3. 


Example 6.3: 


Assume the same setting as introduced in Example 6.1 except that there are neither a mapped 
state signal nor any information about the switches. 

The first iteration of the j-loop of the identification and segmentation algorithm (Fig. 6.8) is 
depicted in the first subplot of Fig. 6.9. 


Both subsystems are considered for verification in each step of the k-loop. It is possible to 
verify subsystem s\“)* fork € {2,3,4,5}. The first result can again be calculated for k = 2 
and the first switch is recognized at k- 2 = 6. This is the first time both dynamic subsystems 
show inconsistency, i.e. Scon,6 = 0. 

Note that the algorithm will break the k-loop at k = 6. However this was not done in Fig. 6.9 
to show that the verification results stay infeasible for all regarded segments in the remaining 
measurement time. 

In the reinitialization step, the detected first segment is deleted from the measurement data. 
This is depicted as shaded area in subfigure 2 and 3 of Fig. 6.9. The second iteration of the 
j-loop verifies subsystem s‘?)* for k € {7,8,9} and detects the next switch at krg = 10. 
Subsystem s\)* is evaluated to be infeasible for the whole available measurement. Again the 
infeasible result for both subsystems are given until the end of the measurement in contrary 
to the genuine break of the k-loop. 

The third segment shows consistency with subsystem s“* fork € {11, 12, 13, 14,15} which 
leads tok, 3 = T. This is feasible with respect to Def. 6.7ff. The subsystem s)* is evaluated 
to be infeasible for segment 3. 

The results can be used to set up the matched state signal 


Q = [1,1,1,1,1,2,2,2,2,1,1,1,1,1,1] (6.66) 


which corresponds again with the ground truth given in (6.40). The successful mapping for 
all time steps implies continuous consistency as there is an unambiguous nominal subsystem 
mapped to each time step. 
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The nominal state machine Z* of Example 6.1 is used here as well. Due to Q = Qmeas; 
all values needed for the verification of the discrete part are the same as in Example 6.1. 
This leads to full discrete consistency between the unsegmented measurement data and the 


specification. 


Finally, full hybrid consistency can also be concluded for this setting. This shows a successful 
verification of a hybrid system only based on interval type measurements of the input-output 
data and the nominal system 57, q 
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Figure 6.9: Verification result for each segment (without segmentation) 
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6.3.1 Convergence of the Identification and Segmentation Algorithm 


To show convergence of the identification and segmentation algorithm it is possible to extract 
more information about the switches from the measurement data by adding a backward iter- 
ation. As a basic property of the identification and segmentation algorithm the true switches 
are always overestimated and never underestimated. This can be shown by considering ex- 
treme values of the additive noise € as follows: 


Noise tending to zero (e — 0) When the noise tends to zero, the interval enclosure nec- 
essary to bound the noise tends to zero, too. Prop. 4.1 has thus to be fulfilled for each 
line of the regressor matrix individually without any deviations. If there is at least 
one measurement point added to the regressor matrix that is generated by a different 
and Prager-Oettli-Distinguishable subsystem, the data is not consistent anymore. The 
switch can thus be detected at the first time step of the new interval, i.e. at the true 
switch with ke; = krj- 


Noise tending to infinity (c — oo) When the noise tends to infinity, the enclosing interval 
width also tends to infinity to contain the noisy data. Therefore the regressor matrix 
provides as well infinite possible entries to fulfill Prop. 4.1 which leads to consistency 
for any measurement data. Hence it is ensured that for the detected switch holds kr, ea 


k 


rj: 


Note that c — oo will also lead to the violation of Prager-Oettli-Segmentability given in 
Assumption 6.4, as well as the violation of the full rank Assumption 3.1. In practice it is thus 
important to ensure a suitable € such that all assumptions of distinguishability, segmentability 
and the rank are met. 

It is not possible to determine the precise amount of time kover = krj — kr,j the switch 
will be overestimated. Nevertheless there will be a gradient from instantaneous detection 
kover = 0 in the case of € — 0 and the maximum overestimation in the case of € — oo. 


To determine the switch more precisely, the detection algorithm can be extended with an an- 
tichronological iteration. The end time of the last segment is thereby given as the time of the 
last measurement kr’ nuien = T by definition. The regressor matrix is set up such that it in- 
cludes the values from a variable time kstart up to the final time T. If the dynamic consistency 


of Prop. 4.1 is fulfilled for all measurement values [m tg 4 (Yanak k tare |? the 


segment is extended backwards in time by reducing kstart <— Kstart — 1. If the segment 
extends over a switch, i.e. kstart < kr,j, data from two distinct subsystems is included in the 
regressor matrix. The detection of this switch is dependent on the same assumptions as in 
forward direction, namely Prager-Oettli-Distinguishability and Prager-Oettli-Segmentability. 
The switches detected backwards in time show the same behavior of overestimation as the 
switches detected in the forward iteration. 


The initial switch is guaranteed to be correct, as it is given by the first measurement data point 
by definition. The results of the backward and forward iteration will occur chronologically 
alternating, see Fig. 6.10. The time between the switch detection in backward direction and 
the switch detection in forward direction is guaranteed to frame the true switching time. 
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The true switching time is thus overapproximated as well in forward direction kr. j,f as in 
backward direction k+ ; ;, leading to a switch segment 


kr € [Kr ke 5,7]. (6.67) 


The measurement data between the switch segments, i.e. for the time steps 


ke [hry f, ke j+1,b] (6.68) 


belong to a so called trust segment. The trust segment is guaranteed to contain only mea- 
surement data generated by a single subsystem. Verification results calculated based on trust 
segments are not disturbed by measurement data generated by other dynamic systems. The 
results of the algorithm converge to the true values for a suitable setting of the enclosed noise 
c and the enclosing interval width ô. 


Switch segment 1 Switch segment 2 
K—— K—————3 K———3 
Trust segment 1 Trust segment 2 Trust segment 3 


Figure 6.10: Example of the alternating occurrence of trust and switch segments 


6.4 Conclusion 


The concept of verification of dynamic systems based on Kaucher interval arithmetic intro- 
duced in Chapter 4 was extended to hybrid systems in this chapter. Therefore a hybrid model 
consisting of a discrete time, value continuous, dynamic system and a discrete event state ma- 
chine were formally introduced. The behavior of the dynamic part is given by parameters 
that are determined by the discrete event states. As long as there is only one discrete state 
active, the dynamic verification can be done as in the nonhybrid case. If the active state 
changes, the procedure has to be adapted. Each segment of measurement data generated by 
a single active subsystem can be verified individually. 


The problem becomes more complex, if it is necessary to determine the segments and the 
active subsystems within the segments. Three different approaches to solve this problem 
were introduced. 


In a first step, the state signal was measured and the given states were assumed to be mapped 
with the specification. In a second step, only the switching times were known and it was nec- 
essary to determine the respective active states. Conditions of Prager-Oettli- 
Distinguishability on the specification and the mappability of the measurement data were 
introduced that ensure the theoretical possibility of determining this information. 
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In a third step, there was no information at all about the switches and the active states. 
Therefore an algorithm was developed that is able to determine both, switches and states, in 
an iterative identification and segmentation procedure. If the introduced property of Prager- 
Oettli-Segmentability holds, this procedure can be applied successfully. 


The verification of the discrete event part can be done in the same way for all three settings. 
The trace given by the mapped state signal is used to set up the measured state machine 
which can be compared easily with the defined state machine. Therefore the set of states and 
the transition function are compared with the nominal values given in the specification. 


This approach was partially developed in [Hen15] and published in [Sch17a]. Also it was ap- 
plied to the practical example of a battery management system [Lem15] and a hybrid braking 
system [Glü17] and it was presented in [Sch16] [Sch18a][Sch19]. 


7 Extended Kaucher Based Guaranteed 
Verification 


This chapter introduces several extensions to the Kaucher based verification framework de- 
veloped in the previous chapters. Therefore the point of view is changed from finding “at 
least one” consistent parameter to calculating the “largest possible set” of parameters within 
the united solution set. This is beneficial as the calculation of the exact solution set is an N P- 
Hard problem as stated in Section 3.1.3. This chapter introduces a method that combines the 
Kaucher based method with optimization techniques to calculate the maximum inner approx- 
imation of the feasible set. 

This approximation is done using different shapes given by the combination of the objective 
function and the constraints. Three different shapes are presented and analyzed. The basic 
approach uses orthogonal enclosure leading to classic interval type borders. 

In a second step, this approach is extended to zonotopic sets to exploit specific properties 
of the measurement data. It is widely known in the interval community that using zono- 
topic sets avoid large underapproximations that arise from axis parallel orthogonal enclosure 
[Jau01][Asa06][Pui06][A1t08 ][Ing09] [Mai16][Roe16][Wan17]. Finally all restrictions on the 
shape of the resulting set are dropped and a general polytope is accepted as result. Although 
this is the most general shape and yields the largest inner enclosures, it is impossible to de- 
scribe the result in terms of interval values, thus limiting the usability of this approach. 

The last section of the chapter introduces an alternative application of the Kaucher based 
verification method in an online setting. Therefore the verification procedure is extended 
to allow an iteratively increasing measurement set. It is thus necessary to provide repeated 
approximations of the solution set as well as repeated verdicts to evaluate the consistency. 
This resembles the situation of an ongoing diagnosis process of a running application. 
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7.1 Solution Set Approximations 


Instead of solving the problem whether a particular (nominal) parameter vector is part of 
the united solution set, the question is now to find the (whole) feasible set given by the 
measurement data. Therefore the constraints cm of the feasibility problem given by the 
measurement data - as defined in Def. 5.6 - are regarded. The feasible set is based only on 
the measurement data and thus no information about the consistency of the VO is included. 
The feasible set is defined as: 


Definition 7.1 (Feasible Set) 
The feasible set F is given by 


f= {o eRratre . c (0) € 0, Vi € (1,2,...,2(T — max (nz, n2) , (7.1) 


Le. the set of all parameters O that fulfill the constraints given by the measurement data as 
defined in Def: 5.6. 


Note that the genuine feasible set here is given according to Def. 5.6 which represents the 
united solution set 5 44. The feasible set is thus not necessarily connected and not necessar- 
ily constrained by borders parallel to the axis (see Section 3.1.2). 


The feasible set can also be defined by a set of vertexes: 


Definition 7.2 (Vertex Based Feasible Set) 


The feasible set F is constrained by the convex hull (see [Bro08, p. 663]), given by the set of 
vertexes V = (V9, Vi, ..., Viv. ik: 


IvI-1 Ivi-1 
F(ysi0em"*:02 V aV Vira: 20^ >) a=1>. (72) 
i=0 i=0 


With this definition it is possible to transfer the setting into an optimization problem based 
on the vertexes of the convex hull. The shape of the resulting approximation of the feasible 
set F becomes a design parameter that is reflected by the objective function and the con- 
straints. The solution of the optimization is thus not necessarily of the same shape as the 
genuine feasible set. 

The approximation of the feasible set can be done in different ways. This thesis uses three dif- 
ferent approaches: hyperrectangular approximation, zonotopic approximation and polytopic 
approximation. All three approaches are introduced in the next sections. 


1 The constraints of the nominal set cy defined in Def. 5.5 are not necessary for the calculation of the feasible 
set. 
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7.1.1 Hyperrectangular Solution Set Approximation 


In the hyperrectangular case, the objective function is set up such that the optimization yields 
the largest hyperrectangle area within the united solution set. This shape represents an in- 
terval type result, constrained parallel to the coordinate axis. 

The solution set is determined using an optimization setting as given in Prop. 7.1. 


Proposition 7.1 (Optimization Based Hyperrectangular Solution Set) 


The hyperrectangular approximation of the feasible set FE is given by the set of vertexes 
VE, The set of vertexes V" is calculated based on the measurement data [U meas, Y meas] 
given for the discrete sampling points k € {1,2,...,T} ofa VO. It is defined as the solution 
ofthe optimization problem 


V= = argmax(J~ (0)) (7.3) 
e 
with the objective function 
nx +n* l 
J^(6)- [] (a - a9) (7.4) 
i=1 
that is subject to the constraints 
c? (e) := 9 0 «o (7.5) 
cm (@) <0 (7.6) 
fori = {1,2,...,n% -- nz). Thereby cp constrains the solution to be proper and cm is 


given by the measurement data according to Def: 5.6. 


Proof: 

The result of the optimization problem (7.3) is given by the set of vertexes Y" that defines 
the hyperrectangular approximation F}. The hyperrectangular approximation is proper in 
all dimensions due to (7.5). The united solution set is given by the measurement data and 
consists of all parameters that fulfill (7.6) according to Def. 5.6. As the optimization problem 
(7.3) is constrained by (7.6), all elements of the resulting set of vertexes are part of the united 
solution set 5 "44 and thus 


FO =F(V") CF Na: (7.7) 
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In general, an n-dimensional parameter space leads to |V7| = 2" different vertexes that 
determine the solution set. In the hyperrectangular case, the set of vertexes V is directly 
given by the interval type parameter vector as illustrated in Example 7.1. 


Example 7.1: 


For ann = 2 dimensional interval type parameter vector 


T 
e- [e^ eo (7.8) 
the resulting rectangle is given by the set of |V| = 2" = 4 vertexes: 
ve - (feto, ao]. fa, so], fo, a] [om s Y. ry 


Without the constraint set cp given in (7.5), the setting can lead to improper solutions. These 
solutions cannot be interpreted and therefore are useless with respect to a real system. This is 
due to the fact that there is no width defined for improper intervals. Nevertheless, the objec- 
tive function is defined on infimum and supremum of the intervals and can thus be evaluated 
even for improper intervals. In case of an even number of parameters being improper, the 
area multiplication in (7.4) yields a positive value. This value can be increased to infinity 
in the "improper direction" leading to impossible results. Such improper solutions are not 
suitable for the verification setting, as the measurement was obtained from a real system, 
with real generating parameters. An illustration of the rectangular inner enclosure is given 
in Example 7.2. 


Example 7.2: 


Assume the united solution of the measurement data to be given as the blue shape in Fig. 7.1. 
The solution is in general not unique because there might be other possible solutions of the 
same size. Two possible area maximal inner enclosures are given by the green rectangles. 
Note that this example is showing the basic concept of the enclosure and thus the denoted 
values are chosen arbitrary. 
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Figure 7.1: Area maximal axis parallel hyperrectangular inner approximation of the united solution set 
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7.1.2 Zonotopic Solution Set Approximation 


In this approach, the united solution set is constrained by a set of hyperstripes in the pa- 
rameter space, generated by the measurement data. The orientation of these hyperstripes 
is similar for consistent measurement data. This leads to a united solution set showing a 
zonotopic shape which is depicted as the shaded area in Fig. 7.2. 


The formal description of the hyperstripes is done based on [Ble11] such that each line of the 
regressor matrix defines a set of feasible parameters F®). Thereby F‘*) denotes the feasible 
set at the specific time step k. The feasible hyperstripes are given by 


FY) = {9 eR" : 5° < Bmeas,k — Ameas,kO X 5°} (7.10) 


with 5° being the width of the interval enclosure, i.e. the maximum absolute sensor fault. 
The hyperstripe can be written in normalized form as 


* * Beas Ameas 
FR) _ fo c gni. | k 5 te x i}, (7.11) 


0.15 — P" 


- ma — a 


0.05 | = 


Figure 7.2: Constraints showing the general zonotopic shape (shaded area) of the united solution set in the case of 
two parameters 


The feasible set F;, that includes all measurement information up to time step k, can be 
determined recursively by intersecting the hyperstripes FR): 


Ra Fane. (7.12) 


Therefore the resulting feasible set for all available measurement data is given by Fr. The 
approximation of this set can be done by a zonotopic set with a very good fit. Definition 7.3 
describes such a general zonotopic set Z as given in [Ble11]. 
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Definition 7.3 (Zonotopic Set) 


A zonotopic set Z is constrained by the convex hull of the set of vertexes 
y= (Vo, Vi,..., Vive a }. The calculation of the vertexes is done by 


y —pP^oHKV = {P+ H"z:zE KV) (7.13) 
with the center of the zonotope P? € R"atrcX1), the radius matrix H € R(atnexV) 


and a unitary box KV composed of an arbitrary number of V unitary interval vectors 


K =|-1, 1j. 


The ®-operator denotes the Minkowski-Sum [Ber08, p. 291] that is used to calculate the 
vertexes of the zonotope. This means that the center P? is added to the given unitary vectors 
(i.e. corners) that are scaled by the respective radius value. The adaption to the zonotopic 
approximation of the feasible set is done in Def. 7.2, leading to the notation F°. The respective 
optimization problem is set up such that a zonotopic approximation F° of the united solution 


Ša; is achieved. 


Proposition 7.2 (Optimization Based Zonotopic Solution Set) 


The zonotopic approximation of the feasible set F° is given by the set of vertexes V^. The 
set of vertexes V? is calculated based on the measurement data [U meas, Ymeas| given for 
the discrete sampling points k € {1,2,...,T} of a VO. It can be computed based on the 
optimal scaling parameter 


o? = argmax (J? (a)), (7.14) 


a 


the initial center P? and the radius matrix H? given by the outer enclosure of the measure- 
ment data. The vertexes of the zonotope V° are calculated by 


y? = P^ oo? HRY = (P? co*H?z:z e KV). (7.15) 
The optimization problem consists of the objective function 
JT? (a) =a (7.16) 
that is subject to the constraints 
em (V°) € 0 (7.17) 


given by the measurement data according to Def. 5.6. 
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Proof: 

The result of the optimization problem (7.14) is given by the set of vertexes V° that defines 
the zonotopic approximation F°. The zonotopic approximation is proper in all dimensions 
as it is based on an initial outer approximation calculated using classic interval arithmetic 
and subsequently scaled using a point real parameter. The united solution set is given by the 
measurement data and consists of all parameters that fulfill (7.17) according to Def. 5.6. As 
the optimization problem (7.14) is constrained by (7.17), all elements of the resulting set of 
vertexes are part of the united solution set 5 ^44 and thus 


Fe =F(V’) CF = Va: (7.18) 


This method has its roots in [Ble11] and [Sch17c] and was developed in [Sch18b]. The initial 
parameters of the optimization are given by the center P? and the radius matrix H? which 
are determined by calculating the outer enclosure as in [Ble11], given in (7.10)-(7.12) extended 
to Kaucher arithmetic notation. 

The initial zonotope to calculate this outer enclosure needs to be chosen suitably. One pos- 
sibility is to use the nominal region O* = [O* — O4, O% + OA] expressed as P" = O* and 
H? = [O5 Itis also possible to calculate the point real central solution ©, using the center 
matrizes Amcas,.c and Byeas,c- The initial zonotope is then given by P? = ©, and H? = Ie 
with an arbitrary small value € > 0. 


Each measurement interval is iteratively interpreted as a hyperstripe containing the possible 
parameters. The intersection between the hyperstripe and the zonotope is calculated and the 
common region is used to calculate the new radius matrix. This procedure leads to a zono- 
topic outer enclosure of the feasible parameter set. Due to the repeated calculation of outer 
enclosures, the area of the zonotope might grow with the considered measurement data. 
Starting from the final outer enclosure that frames the intersection of all available measure- 
ment data, the scaling factor o is minimized until all vertexes of the zonotope are part of the 
united solution. This can be checked using Prop. 4.1. 

The center point of the zonotope is not moved by the shrinking procedure. The idea of maxi- 
mum possible parameter variability within the zonotope is realized by the objective function 
(7.16) that maximizes the scaling factor a. 

Additional assumptions on the optimization problem are: 


+ The interval solution has to be bounded to one orthant. 


* All values of the input signal need to have the same sign, either all positive or all 
negative. 


The first assumption is due to the fact that intervals containing zero can be interpreted erro- 
neously as inverse elements and thus cancel the influence of some parameters. The second 
assumption is necessary to prevent increasing intervals that may arise even when Kaucher 
interval arithmetic is used. 

A sketch of a zonotopic inner approximation is given in Example 7.3. 
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Example 7.3: 


Assume the united solution of the measurement data to be given as the blue shape in Fig. 7.3. 
The center of the zonotope is depicted by the green circle. It is given by the center of the outer 
enclosure that is used as initial zonotope. The algorithm does not change the center, leading 


to the area 


maximal inner enclosure given by the green zonotope. 


Note that there might be solutions of equal or larger size possible for a different zonotope 
center. However the choice of an optimal center is not included in the current version of the 


algorithm. 


0.20 


The denoted values are chosen arbitrary in this example. 
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Figure 7.3: Exemplary area maximal zonotopic approximation of the united solution set 


7.1.3 Polytopic Solution Set Approximation 


The most general enclosure is given by the shape of a convex polytope. However in this case 
the solution is represented by a list of points instead of a mathematical description such as 
intervals (in the hyperrectangular case) or center and radius matrix (in the zonotopic case). 
The only conditions on a valid polytopic enclosure consists of the requirement that points of 
the list are part of the united solution set 5 "44 according to Prop. 4.1 and that the resulting 
polytope is convex. The computational effort of the polytopic approximation is in the same 
order of magnitude as the computational effort of the hyperrectangular approximation. The 
respective solution of the optimization problem is given in proposition Prop. 7.3. 
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Proposition 7.3 (Optimization Based Polytopic Solution Set) 

The polytopic approximation of the feasible set F? is given by the set of vertexes V°. The 
set of vertexes V? is calculated based on the measurement data [U meas, Ymeas| given for 
the discrete sampling points k € {1,2,...,T} ofa VO. 

It is defined as the solution of the optimization problem 


y? = argmax (J? (©°)) (7.19) 
o9 
with the objective function 
J? (O°) = area(0?). (7.20) 


The function area (-) calculates the hypervolume of a polytope given by a list of points O?. 
All points of the list, i.e. vertexes of the set F°, are subject to the constraints 


em (O°) € 0 (7.21) 


with cm given by the measurement data according to Def. 5.6. 


Proof: 

The result of the optimization problem (7.19) is given by the set of vertexes Y? that defines 
the polytopic approximation F°. The polytopic approximation is proper in all dimensions 
as it is the convex hull of V°. The united solution set is given by the measurement data and 
consists of all parameters that fulfill (7.21) according to Def. 5.6. As the optimization prob- 
lem (7.19) is constrained by (7.21), all elements of the resulting set of vertexes are part of the 
united solution set 5 ^44 and thus 


F?LJF(V)cF-u. (7.22) 


The choice of the vertexes of the initial list is of minor importance. Possible choices are the 
vertexes of the nominal set, the vertexes of a zonotopic outer enclosure or the vertexes given 
by the central solution disturbed by a small parameter e > 0. 
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Example 7.4: 

Assume the united solution of the measurement data to be given as the blue shape in Fig. 7.4. 
A possible inner enclosure is given by the green polytope. 

Note that this solution is not unique. There might be other possible solutions of the same size. 
Again, this example is showing the basic concept of the enclosure and thus the denoted values 
are chosen arbitrary. 
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Figure 7.4: Exemplary polytopic inner approximation of the united solution set 


7.2 Kaucher Based Diagnosis 


The consistency for interval type systems was introduced in Chapter 5. Prop. 5.4 focuses on 
basic consistency and includes the main result that forms the foundation of the considera- 
tions in this chapter. All other assumptions and definitions are assumed to be fulfilled and 
applicable as well. 

In case the calculation of the verification method can be fast enough with respect to the re- 
garded dynamic system, the method can be applied in an online setting to tackle the diagno- 
sis problem. In the resulting diagnosis setting, the interpretation of the problem includes the 
temporary feasible set F+ that is approximated using one of the three approximation shapes 
introduced in Section 7.1. The intersection between the temporary feasible set Fj, and the 
nominal set N forms an inner enclosure of the consistent set for the whole measurement 
time according to Sec. 5.2. 
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Proposition 7.4 (Approximation Based Basic Consistency) 


The input-output behavior given by all available interval type enclosures of T measurement 
values 


T T 
[U eas; Y meas] = er a en pci (7.23) 
is basic consistent with the nominal system specified by an interval type specification 
Sj = {0*, ni no Vanity You), (7.24) 
if the intersection between nominal set N and temporary feasible set F;, is nonempty, i.e. 


NOF, EI. (7.25) 


Proof: 

The nominal set N consists of all parameters within the specification ©*. Basic consistency 
according to Prop. 5.4 is given, if there is at least one parameter © € ©* within the specified 
parameter set that is able to explain the measurement data. 

According to Def. 7.1, all parameters of the feasible set Fẹ are able to explain the measure- 
ment data for k € (1,2,..., T}. 

The intersection between the nominal set N and the feasible set Fp contains parameters that 
are both, part of the nominal set and part of the feasible set. If this intersection is nonempty, 


there is at least one parameter that is consistent according to Prop. 4.1. 


It is in general not possible that the method yields full consistency in the diagnosis setting. 
However, the approximation methods described in the previous section can be used to verify 
dynamic systems given by an interval type specification 57.4 

In the diagnosis setting, an approximation of the feasible set is used to calculate an inner en- 
closure of the consistent set. The verification result is guaranteed in the sense of this thesis as 
it is based on Kaucher arithmetic. The definition of consistency uses the united solution set 
as given in Def. 3.1. The choice of vertex points used during the verification is done based on 
the introduced optimization procedure hence ensuring an effective coverage of the feasible 
area given by the measurement data. The shape of the calculated solution is determined di- 
rectly by the setup of the objective function and the additional constraints. Furthermore, the 
calculated inner approximation aims at finding an area maximal representation of the united 
solution set in terms of the defined approximation shape. Thus the feasible set represents the 
maximum possible parameter variability in the given setting that is guaranteed to be able to 
explain the measurement data. 


14 This method is also applicable to a point real specification. Nevertheless a point real specification can be verified 


directly by Prop. 4.1 and does not need the optimization based extensions in order to approximate the feasible 
set. 
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It is assumed that the diagnosis algorithm is running in parallel with the VO. The diagnosis 
system is supplied with a new measurement value in each sampling cycle which is used to 
calculate a temporary result. This leads to the diagnosis algorithm as given in Fig. 7.5. 


Initialization of 
the Iterative 
Algorithm 


Get Nominal 
Set N from 


Specification 


Calculate Ap- 
proximation of 
Feasible Set JF; 


Intersect Approxi- 
mation of Feasible 


Set 7, with 
Nominal Set N 


Measurement 


Measurement is . : . 
: Y no s might be in- 
basic consistent 


with Specification 


consistent with 
Specification 


Figure 7.5: Diagnosis algorithm 


A VO that is basic consistent for all measurement data at the final time T, considering all 
measurement data with k < T according to Prop. 7.4, is assumed to be also basic consistent 
fort <T. 

If the algorithm yields basic consistency between the measurement data and the specifica- 
tion, the verdict is guaranteed to be correct. Therefore there are no hidden faults i.e. no type 
I errors possible. 

If the algorithm yields inconsistency, the verdict might be erroneous, i.e. showing a false 
alarm. This is due to the used inner approximation of the feasible set which does not neces- 
sarily cover all feasible parameters. 


Three exemplary settings and the resulting solution sets are depicted in Fig. 7.6. The black 
lines depict the constraints given by the measurement data that frame the united solution. 
The nominal parameters are given by the black dashed rectangle. 
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Figure 7.6: Exemplary results for the three different approximation shapes 


The first plot shows an inner approximation using an rectangular shape. It can be seen that 
the location and shape of the rectangle is not unique and that other rectangles of the same size 
might be possible within the united solution. The second picture shows the approximation 
by a zonotope. The general fit of the zonotope is better as it can follow the contour of the 
constraints to some extent. The third figure shows the approximation by a polytope. This 
setting shows the best fit, even though the resulting shape is given by a list of four points 
instead of an algebraic description. All three examples show an intersection between the 
nominal set and the approximation of the united solution set, given by the shaded areas. The 
regarded example setting thus shows a situation with guaranteed basic consistency between 
measurement and specification. 
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7.2.1 The Center Misplacement Effect 


The zonotopic approximation leads to the best tradeoff between accuracy and ease of de- 
scription. This shape is chosen for the further explanations of the diagnosis approach. As 
already mentioned, false alarms are possible for all regarded shapes due to the usage of inner 
approximation of the feasible set. 

In the special case of a zonotopic approximation, false alarms can also result from the so 
called Center MisPlacement Effect (CMP). CMP denotes the effect of disregarding feasible 
results due to a bad center and shape of the zonotope. 

An exemplary situation showing a CMP effect is given in Fig. 7.7. 


0.050 


Figure 7.7: Center misplacement effect 


The blue lines depict the constraints given by the measurement data. The center of the zono- 
tope is given by the green circle within the green zonotope. Center and shape of the zonotope 
were calculated from an outer approximation of the solution set as in [Ble11]. Afterwards 
the scaling parameter o was adapted such that all vertexes of the zonotope are within the 
united solution of the problem, given by the shaded area. Temporary feasible vertexes of the 
zonotope that were checked during this optimization procedure are given by green crosses." 
Due to the location of the center and the shape of the zonotope, the scaling factor a needs to 
be very small to ensure that the vertexes V2 and V3 stay within the united solution set. 

The resulting zonotope does not intersect with the nominal set, from which the lower right 
corner is depicted by the dashed black square in the left of the figure. 

Even though the final vertex is rejected there are consistent intermediate vertexes (green 
crosses within the nominal set) that are rejected because the respective temporary zonotope 
was rejected as not all vertexes were part of the united solution for this value of the scaling 
parameter a. 


Certain measures can be taken to avoid CMP. The most straight forward is to move the center 
of the zonotope or to change its shape. However, in general it is not trivial to chose a suitable 
change in center and shape algorithmically. 


D Temporary vertexes of the zonotope that were checked during the optimization procedure and found to be 


infeasible are marked as red circles. 
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As the method is also designed for higher dimensions and should work as automatically 
and autonomously as possible, it is neither possible nor suitable to visualize the setting and 
request a human operator to adapt the center or the shape. 

An algorithmic workaround is given by the use of intermediate vertexes that are checked 
during the optimization. The optimization can be stopped as soon as there is at least one 
feasible intermediate vertex detected, i.e. if a green cross is evaluated. Nevertheless this will 
lead to point-wise results instead of a feasible set of a given shape. An effective way to avoid 
CMP is provided by the enclosure of the nominal set in the constraints of the optimization 
problem. However, this will also pose higher restrictions on the feasible set. 


7.3 Conclusion 


This chapter introduced extensions to the Kaucher based verification method developed in 
the previous chapters. 

The first extension is given by calculating the largest inner approximation of the united so- 
lution set. This was done using an optimization setting. The precise setup of the objec- 
tive function and the constraints determines the resulting shape of the so called feasible set. 
The feasible set was approximated using three different geometric shapes (hyperrectangle, 
zonotope, polytope). The resulting set is an inner approximation of the united solution set, 
meaning that there are several equivalent approximations within the united solution set. 


The second extension is given by the application of the verification method in a diagnosis 
setting. The diagnosis algorithm was introduced in an iterative setting, calculating temporary 
results for each sampling cycle. Due to the inner approximation used, the results can show a 
false alarm. A vivid cause of a false alarm is given by the center misplacement effect, that can 
occur in case of a zonotopic approximation. The applicability of the Kaucher based diagnosis 
method depends on the specific settings of the regarded system. Therefore it is necessary to 
evaluate each system individually before applying Kaucher based diagnosis. 


8 Application and Results 


This chapter presents the application of the developed methods to tank systems with a vary- 
ing number of tanks and an adjustable set of connections. The settings are analyzed in sim- 
ulation and practice. 

Tank systems form a class of wide spread theoretic control applications. The basic setting is 
usually given by one or more tanks with a nominal outflow and a controllable inflow. The goal 
is to adjust or maintain a nominal height of the fluid in a tank. The process can be disturbed 
by additional leakages or inflows or by congestions in the in- or outflow pipes. Depending 
on the specific setup, cross flows between tanks are possible. Those cross flows are in gen- 
eral dependent on the current filling level of the concerned tanks. All measurement data is 
obtained using a real three-tank process available at the Institute of Control Systems (IRS) at 
the Karlsruhe Institute of Technology (KIT). A picture of the lab setting is shown in Fig. 8.1. 
The algorithmic calculations are done on a Lenovo ThinkPad T460s powered by 
an Intel® Core™ i7-6600 CPU using 12GB main memory. The implementation was 
done in Matlab® 2012b. 


Figure 8.1: Three-tank process lab setting at the Institute of Control Systems (IRS) 
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The following subsections set up the dynamic models and introduce the properties and pos- 
sibilities of the different scenarios. First the methods of Chapter 5 “Guaranteed Verification 
of Interval Type Systems” are applied to real measurement data obtained from a single-tank 
process. 

Second, the application is extended to a two-tank setting showing hybrid behavior. A mapped 
state signal is used to demonstrate the basic functionality of the method developed in Chap- 
ter 6 “Guaranteed Verification of Hybrid Systems”. The results are discussed and interpreted. 
Third, the diagnosis method developed in Chapter 7 “Extended Kaucher Based Guaranteed 
Verification” is applied to simulation data of a four-tank process. Three different fault types 
are used to demonstrate the fault detection properties of the Kaucher based method. Several 
fault intensities demonstrate the possibility of the method to detect even very small faults. 
Finally the diagnosis method is applied to real measurement data provided by a single-tank 
process. It is shown that all of the regarded faults can be detected successfully using the 
introduced methods. 


8.1 Application: Guaranteed Verification for Interval 
Type Systems (Single-Tank) 


The basic setting is given by a single-tank process which is sketched in Fig. 8.2. 


Pump 2 va 


Nominal outflow valve Vout2 | 


Nominal outflow floz 


Figure 8.2: Sketch of the single-tank 


It is possible to set up a single-tank system consisting only of tank 2 by closing the respective 
valves in the three-tank lab system.!° All valves in the system are binary valves which are 
only open or closed. The height ha of tank 2 is measured, as well as the flow v2 of pump 2. 


16 Note that the number of the tanks in the lab setting is (from left to right) 1 - 3 - 2. 
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The outflow of each tank is governed by the formula of Torricelli [Tip00, p. 360] which leads 
to the nonlinear time continuous dynamic of a single-tank 


dha(t 1 1 
P e. =, av 2gha(t) +7 - Yava(t). (8.1) 
outflow inflow pump 


with the outflow pipe cross section a», the tank cross section A», the gravitational force g 
and the constant Ya according to Appendix G, Tab. G.1. The following simplifications and 
the resulting model equations are based on the considerations in [Ble11]. 


The model is discretized using the Euler method with sampling time At leading to 


a 
hə k =hop—1 — P 2gh» y 1 At + Atari AL + ea, (8.2) 
2 2 


where e2 ẹ is the additive disturbance including sensor and discretization faults. Equation 
(8.2) is reformulated to the pseudo linear regressor form 


Pkk = Yk (8.3) 
with 
Pk = ha,r-ı (8.4) 
m=i EN p o. m 
Yk = ha — AL vzet. (8.6) 


It can be seen that the parameter 0, given in (8.5) is depending on the height A» ;..; which 
renders it time variant. The range of the time variant parameter can be interpreted as an 
interval set that includes all possible parameter values as well as some spurious solutions. 
The interval enclosure of the parameter is given by the bounding box of the time variant 
parameter. 


It is possible to calculate the interval enclosure of the time variant parameter 6; for a spe- 
cific nominal setting (ea, = 0). This setting consists of a given operation range ha € 
[ho min; R2,maz] and a fixed sampling time At. The calculated parameter interval can then 
be used as the nominal set in further considerations. 


Based on the tank properties given in Appendix G, Tab. G.1 it is possible to set up the pa- 
rameter range for a nominal and a faulty tank configuration. To realize the faulty behavior 
all available valves are opened. This means the leakage outflow valve, the lower connec- 
tion valves and the nominal outflow valves. The resulting height dependent parameter 0x is 
depicted in Fig. 8.3 and can be used a priori to reason about detectability. 
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Height hy in cm 


— Nominal — Faulty 


Figure 8.3: Values of 0; depending on ha for different outflow configurations 


It can be seen that there is a gap between the nominal behavior (green line) and the faulty 
behavior using all possible outflows (red line). This leads to the hypothesis that it is possible 
to separate the nominal behavior from the faulty tank configuration using the method from 
Chapter 5. Due to the upper valve position hyu = 30cm, the depicted values are valid in the 
case of a closed upper connection valve v3,, only. The nominal parameter 0* = [0.93, 0.98] 
is used for tank levels in the range h2 € [5, 30] cm. 


The hypothesis is verified using the following steps: First, a nominal description i.e. a specifi- 
cation, of the system in the necessary ARX form is set up. Second, the system is implemented 
ie. the real tank is manufactured and used to collect measurement data. Third, a faulty ver- 
sion of the system is implemented, i.e. unspecified outflows are added to the tank by opening 
the respective valves. 

The collected measurement data for ha is then enclosed by interval values 


ha = [h2 (1 — 67,), ha (1 + õn, )] (8.7) 


and used to verify the correct system with respect to the nominal parameters. Finally, mea- 
surement data from the faulty system is used to show that it is not possible to verify the faulty 
behavior using the nominal parameters. 


8.1 Application: Guaranteed Verification for Interval Type Systems (Single-Tank) 113 


The initial height is set to h2 min = 5cm, the nominal outflow valve Vout2 is open and pump 2 
is constantly running, providing the maximal inflow of v = 6.51/m. The resulting nominal 
setting is depicted in Fig. 8.4. It can be seen that it is possible to verify the measurement of the 
fault free system supposing a relative fault of 6; = 0.02 for the measurement signal of the 
height ha. The basic consistency as introduced in Chapter 5 is used to calculate the results. 
The calculation time for the algorithm is Teate = 16.8s which is less than the duration of the 
experiment T' = 60s. 
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Figure 8.4: Implementation of the single-tank setting that is basic consistent with the specification 
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Now a faulty implementation of the tank system is considered. Therefore the faulty setting 
is realized by opening the leakage outflow, the lower connection and the nominal outflow 
valves which leads to a major change in the system dynamics. 

The changed dynamics are not able to reach the nominal final height of hz 69 = 29.1cm as the 
maximum pump flow cannot compensate the additional outflow. The resulting measurement 
data is depicted in Fig. 8.5. It is again analyzed using a relative fault of ô, = 0.02. Fig. 8.5 
shows that it is not possible to verify the faulty implementation for k > 4. The necessary 
calculation time is Teale = 47.58. 
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Figure 8.5: Implementation of the single-tank setting that is inconsistent with the specification due to additional 
outflows 
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The change in the system dynamics has to be strong to prevent the verification of the faulty 
system. If there is only a slight change in the system dynamics, the faulty system is verified 
because the new behavior can be explained with the range of the nominal parameter set. This 
is the case, if a faulty system is implemented consisting of other combinations of outflows 
and connection valves except the introduced setting. 

When having a closer look on the tank parameters, this is rather intuitive. The connection 
valves cross section a32, and a32,, are the same as the nominal outflow aa. The leakage aleak 
is only slightly bigger that the nominal outflow. Thus each individual valve leads to none, 
respectively very slight variations with respect to the nominal behavior. Any combination of 
two outflows is also verified using the original setting. 


This concludes the first application example concerning the guaranteed verification for in- 
terval type systems as introduced in Chapter 5. It is clear that the results are depending on 
the value of the used relative fault ô; . Increasing the fault increases the variety of enclosed 
trajectories and can thus be interpreted as increasing the available system behavior. This 
leads to a higher chance to achieve nominal behavior enclosed in the measurement data and 


thus to verify the system. 


8.2 Application: Guaranteed Verification for Hybrid 
Systems (Two-Tank) 


The single-tank-system can be extended by adding another tank, connected via two horizon- 
tal valves. There is no additional pump and thus no external inflow to the new tank. However 
the new tank has a nominal outflow and it is possible to open and close the connection valves. 
A schematic description is depicted in Fig. 8.6. 

The dynamic of the main tank needs to be extended with the flows induced by the new tank. 
Those flows depend on the height differences between the two levels, related to the static 
height of the lower and upper valves A32; and h32u: 


A534 = max(ha k, haoj)) — max(h3 k, haaj) (8.8) 


Ao 3,u,k = max(ha k, hasu) — max(h3 k, hao). (8.9) 


The resulting model is again discretized using the Euler method with sampling time At 


1 L 
han = hai — As az y/ 2gh2,k—1At -77 sign(A» 5.1.4 1)a3214/ 2g|^2,3,1,k—1|At 
= re 


2X. 
outflow lower cross flow tank 3 


T. 1 
EET sign(A2,3,u,k—1)432u4/ 2g|A2,3,u,k-1|At T ¥2V2,~-1At +€2,% 
As Ag CLR o 


————————————— 


upper cross flow tank 3 


inflow pump 2 


(8.10) 


with all values according to Appendix G, Tab. G.1. 
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Figure 8.6: Schematic sketch of the two-tank experiment 


The cross sections of both tanks are the same, i.e. Ag = As, leading to symmetric height 
changes induced by the cross flow. The term e» ;; is the additive disturbance including sensor 
and discretization fault. The system description can again be transformed to the pseudo linear 
regressor form (8.3) with 


Pk = [ha,r-ı, |A2,3,1,0-1|, |Ae,3,u,4—11] (8.11) 
E 
©; = e, 9), 9) (8.12) 
Uk = hap — Fhe At. (8.13) 
2 


The elements of the parameter vector ©; are given by 


a2 2g 
gU =1 Jd At+ eh, 8.14 
n A» V har-ı e ( ) 


a 2 
9% = sign(A» ENI k—1) 3al 9 At4 e3 k (8.15) 
= Ag WV |Aa,3,1,x-ıl i 
3 A A32u 2g m 
0i = sign(A» 3.u k—1) At H Cok (8.16) 
= Ag V |Ao,3,u,k—1| i 


with unknown composition of the fault ea, x = €2,«/n2, 1 + ©2,%/|A2,3,1,6-1| + £2, | Asa uihil 
The parameters 92 and ee given in (8.15) and (8.16) show singularities in case |A» 5.1. .1| 
or |A2 3.u,«-ı| approach zero. 
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Therefore the enclosing intervals 0?) and 6‘) become very large when the operation range 
[h2 min, R2,maz] is including or close to the height of the valve h32,,. It is thus necessary to 
chose the operation range with a sufficient distance to 32, to achieve meaningful parameter 
intervals. 


8.2.1 Measurement With Mapped State Signal 


The first approach for the verification of the resulting hybrid system was introduced in Sec- 
tion 6.1 and used a so called “mapped state signal” according to Definition 6.10. 

This means that the exact switching times and the respective active subsystems are known 
correctly. Therefore the hybrid verification task is reduced to sequential verification of the 
subsystems present in the measurement data. 

The general system behavior given in equation (8.10) is therefore transferred to a more spe- 
cific setting. Tank 3 is assumed to be empty (h3,1 = Ocm) with its nominal outflow valve 
open and the lower connection valve v32; closed. The upper connection valve v32, is open, 
as well as the nominal outflow valve of tank 2. 

The resulting hybrid scenario consists of two states: State 1 is active if ho, < hgau, ie. 
the upper connection valve does not influence the system dynamics. State 2 is active for 
ha,r > hs2u, ie. an additional outflow is given through the upper connection valve. 


The hybrid scenario is described as follows: 

In state 1, starting at ha,ı = 5cm, pump 2 is used to fill tank 2. The nominal behavior of tank 
2 is identical with the single-tank behavior described in Section 8.1. 

State 2 is reached, when ha rises above the height of the upper connection valve, i.e. 
har > haa, = 30cm. Thus the system dynamic changes due to the additional cross flow 
from tank 2 to tank 3 through the connection valve. 


The resulting dynamic of state 2 is given by 


1 
har = har-ı — 2, | @v 2gho,n—1At + azau\/ 2glha2,r-ı — hazu| At — avo 1A +e2,% 
PM (a piu MT Exe 


SZ, 


nominal outflow upper cross flow to tank 3 inflow by pump 2 


(8.17) 


In state 2 the level of tank 2 is always higher than the upper connection valve i.e. 
h2 k—1 > haau. Therefore the absolute value operator |: | on (ho,..1 — hg2u) can be dropped 
and the singularity present in (8.16) is avoided. 
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The upper cross flow to tank 3, given in (8.17) can thus be reformulated: 


as2u / 29 (ha s — ha2u) At (8.18) 


2 
= = a32u(ha, k—1 7 hao) ai (8.19) 


= aasu 324, 8.20 
32u12,k—14] PE A 32u/132u4 | 7— — —L.—3 (E A (8.20) 


=h a32u At — azauhaau At 
Tm ( à eo did y hl, m um 


(8.21) 
Finally, the model for state 2 is given in in pseudo linear regressor form (8.3) with 
Pk = ha,k-a1 (8.22) 
a? | 2g — | 
0, —1 At 
: Ag V haa (ha, k-1— a 
Q32u €2,k 
u ae, A 8.23 
A, M Hour Rea Rees (29 
Uk = ha, — Ari At, (8.24) 
2 


The nominal parameter 0j, can again be determined depending on the level h2, ķ—1 for differ- 
ent outflow configurations. The resulting parameter ranges are depicted in Fig. 8.7. and are 
used to determine the nominal values of the system parameters. 
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Height hə in cm 


— Nominal — Faulty 


Figure 8.7: Values of 6; depending on ha for different outflow configurations at state 1 (left) and state 2 (right) 
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State 1 is assigned with the same values as in Section 8.1 i.e. 
0(1)* — [0.93, 0.98]. (8.25) 


In state 2, i.e. starting from h2 ķ—1 = 30cm, the time variant parameter for nominal outflow 
only can be enclosed by the interval 


0(2)* = [0.965, 0.974]. (8.26) 


The nominal parameters are now used to verify the hybrid setting. Again the collected mea- 
surement data of the height ha is enclosed by interval values 


ho k = [h2 (1 — 05,), hau (1 + 65,)] (8.27) 


with ôr, = 0.15. 


The upper part of Fig. 8.8 shows a hybrid test run, starting with tank 3 being empty and the 
level of tank 2 being at hg, = 5cm. The level of tank 2 is increased using pump 2 and thus 
crosses the valve height h32,, at k = 64 seconds. The tank dynamics are changed by the 
crossing as the upper connection valve now acts as an additional outflow of tank 2. 

The verification results are depicted in the lower part of Fig. 8.8. State 1 can be verified as 
long as the level in tank 2 is below the connection valve, i.e. hz, < h32u. Afterwards, state 2 
is verified until the end of the measurement sequence. 


Note that in this case it is not necessary to perform cross validation, i.e. applying the param- 
eters of state 2 to measurement from state 1 and vice versa. This is due to the fact that the 
real switching time and the respective active state are given by the mapped state signal. 
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Figure 8.8: Verification result in the hybrid case for state 1 and state 2 with a mapped state signal 
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8.2.2 Measurement Without Mapped State Signal 


The introduced setting is now generalized by omitting the mapped state signal, still assuming 
known switching times. This means that it is still known there is a state change at k = 64 in 
the scenario, but now it is unknown whether the system switches from state 1 to state 2 or 
vice versa. Therefore all segments have to be analyzed twice, using the nominal parameters 
of both states. This cross verification is used to determine the active state of the respective 
subsystem. 

To ensure successful cross validation, Prager-Oettli-Distinguishability as given in Assump- 
tion 6.2 has to be fulfilled. Therefore the nominal parameters of state 1 and state 2 are in- 
vestigated. It is obvious that the nominal parameters of state 2 are a subset of the nominal 
parameters of state 1, as 


0* (1) = [0.93, 0.98] > [0.965, 0.974] = 0* (2) (8.28) 
which leads to the existence of a common parameter 


Ož m = 0* (1) n 6* (2) = (0.965, 0.974]. (8.29) 


com 


Therefore all parameters that are within the set of common parameters 0%, are consistent 
with both states by definition. 

It is hence not possible to distinguish the states as Prager-Oettli-Distinguishability is not 
fulfilled in the current setting. Since this property is the main preliminary of hybrid verifica- 
tion without mapped state signal it is formally impossible to demonstrate the viability of the 
method using this setting. 


The performance of the algorithm in case of fulfilled Prager-Oettli-Distinguishability was 
demonstrated in Example 6.2. Furthermore it was shown that it is possible to segment and 
verify the measurement data even without information about the switches if Prager-Oettli- 
Segmentability holds. The respective setting and the results are given in Example 6.3. 
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8.3 Simulation: Diagnosis By Kaucher Based Guaranteed 
Verification (Four-Tank) 


Figure 8.9: Schematic view of the used four-tank system 


The regarded setting is now changed to a slightly different four-tank setup. The four-tank 
process is an established benchmark in literature and was proposed by [Joh00]. Here it is 
used to apply the diagnosis method introduced in Chapter 7. This application was published 
and presented in [Sch18b]. 

The four-tank setting is depicted in Fig. 8.9. For symmetry reasons the setting can be reduced 
to tank 1 and tank 3. 

The dynamic of tank 1 is chosen to be the objective of the verification. The heights hı and ha 
of both tanks are measured, as well as the on/off signal v, of pump 1. The flow of pump 1 is 
split by the input valve vj; leading to v;;1 = 0.7 of the flow going to tank 1 and (1 — v;nı) 
of the flow going to tank 4. The flows from pump 1 as well as from tank 3 are considered as 
inputs. The respective equations are similar to (8.2), now taking into account an additional 
inflow depending on hs. 


All simplifications and the resulting model equations are based on the considerations of 
[Ble11]. 
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Discretization using Euler Method and sampling time At leads to 


a1 1 1 
hia =hır-ı — A, dv 2ghi ki At tI az y/2gh3 k—-1At m Uin1^f1U1,k—14M +€1, k 
] ——— 3 N—— 1 —— 
outflow inflow from tank 3 inflow by pump 1 


(8.30) 


with parameters according to Appendix G, Tab. G.2. The additive disturbance eı ; includes 
sensor and discretization faults. 
The pseudo linear regressor form (8.3) is now given by 


Yk = hi — pm Ui, k—1 4A (8.31) 

1 
Pr = [hı,r-ı h3,k-1] (8.32) 

T 
9, = [ec ee (8.33) 

with the time variant parameters 
o =1- ^ 29 At (8.34) 
Ai \ hı,r-ı 

2 

Je 9 At (8.35) 


It is assumed that the operation range of the tank system is hı € [2, 10.5] cm and 
ha. € [1, 15] cm which leads to 9)* = [0.921, 0.965] and 6°?)* = [0.029, 0.112]. 

The resulting midpoint radius expressions ofthe parameters are 6 (?* = 0.943, e — 0.022 
and 6°)* = 0.0705, 9 2* = 0.0415. 

The optimization based diagnosis approach using a zonotopic approximation is chosen in this 
example. Therefore the measurement data and the nominal parameter set are used to set up 
the constraints of the optimization problem. 

The nominal feasible parameter box ©* is used to build the initial zonotope with: 


T 
po | gi» g | (8.36) 
(1) 
H? = ÜA p f (8.37) 
0 o% 


Afterwards the outer enclosure of the intersection between initial zonotope and measurement 
data is calculated using (7.10)-(7.12). The resulting zonotope is used as starting point PP, HQ 
of the optimization problem. The solution of the optimization problem is thus a zonotopic 
approximation of the united solution set given by Z. All parameter vectors included in the 
optimal solution set Z are solutions of the ILES (3.46). If the intersection of the specification 
and measurement is nonempty, the algorithm calculates a feasible set in this area. The results 
and limitations of the approach are demonstrated in the following, using several different 
settings. 
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8.3.1 Fault Free Setting 


First a fault free scenario is given as depicted in Fig. 8.10. The scenario includes parts with 
pump on and off and thus shows a variety of different water level dynamics both in tank 1 
and tank 3. The measurement data of hı and hg are enclosed using intervals with radius 


2 =ô}, = 0.05cm leading to the interval values 
L 3 


hi = [hi — Sk hie + OF, | (8.38) 
hs p = [ha — h, han + 0p]. (8.39) 


The results of the optimization based zonotopic method are given in subplot 4 of Fig. 8.10. 
The algorithm calculated a feasible set of parameters for the time segments [1, kena] with 
kena € [[1, 1297], [1572, 2000]]. However, time segments starting in the beginning and 
ending in kena € (1298, 1571] are not verified. Therefore there is temporarily no consistency 
according to Prop. 7.4 given in this segment. This is due to the CMP effect as introduced in 
Section 7.2.1. 

A detailed view on the relevant time instants is given in Fig. 8.11 which displays the change 
of the result from consistent to inconsistent and back. 


The constraints given by the measurement data are depicted by the blue lines in Fig. 8.11, the 
nominal set by the shaded area and the zonotopic approximation of the united solution set is 
shown as the green zonotope. 

At k — 1297, the measurement data is proven to be basic consistent with the specification as 
there is a nonempty intersection between nominal set and the approximation of the united 
solution set (orange part of the zonotope in Fig. 8.11, subfigure 1). 

At k — 1298, the verification result is inconsistent for the first time. However, there is still 
a feasible region within the nominal parameter set, shown by the green crosses that depict 
vertexes fulfilling Prop. 4.1 (Fig. 8.11, subfigure 2). Those points were used by the algorithm 
while calculating a suitable factor o. However there is no intersection between the final 
green zonotope and the shaded nominal region. This behavior reflects exactly the definition 
of the CMP effect introduced in Section 7.2.1. The CMP effect is observable until k — 1571, 
(Fig. 8.11, subfigure 3). 

Starting from k — 1572, additional constraints given by new measurement data are taken 
into account. Therefore, the center and shape of the outer enclosure is changed. This results 
in a zonotopic approximation of the united solution set that provides a nonempty intersection 
with the nominal set again (Fig. 8.11, subfigure 4), leading to a successful verification. 

As those later results are calculated based on all measurement data, including the possibly 
inconsistent times k € [1298, 1571], the results showing the CMP effect are corrected and 
the verification result for k > 1571 can be generalized for all k € [1, 2000]. 
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Figure 8.10: Verification result for the fault free setting 
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The method is applied iteratively to calculate an individual verification result for each time 
step k. The calculations are done offline and the necessary calculation time for the complete 
data set is Teate = 703.88 < 2000s = T. 

Regarding the calculation time of each measurement time step shows that Teatc,k < 1s which 
means that the method might in general be real time capable. However, the algorithm is based 
on an optimization procedure with non-deterministic runtime. Therefore special measures 
need to be taken to ensure deterministic runtime of each step. For an initial estimate the 
average runtime of the optimization algorithm is used by regarding the total runtime of the 
method throughout this chapter. 


8.3.2 Additive Faults 


There are several phenomena that can lead to an additive fault of an sensor. A possible 
sensor fault is called “freeze”, when the sensor will return a fixed constant value. Another 
common sensor fault is “offset”, which means that the sensor will add a constant bias to the 
true measurement value. A third additive sensor property is the specific sensor noise. Noise 
is not regarded as an effect to be detected here. However it is crucial to know the sensor noise 
precisely to choose the bound of the interval enclosure correctly. 


Given the correct faultless but noisy sensor data 5;, a freeze fault of value f, occurring at 
time kerr, can be expressed as follows: 


sk Whe 1, kapp = 1] 


(8.40) 
f; Yk € [kerr, T]. 


Sf,k = 


The measurement values of hı and hg are enclosed using 6;,, = 6;,, = 0.05cm. The case of 
a freeze fault of fr = 6.0cm on measurement hj at kerr = 650 is depicted in Fig. 8.12. It can 
be seen that the detection time is equal to the fault time kaet = 650 = kerr which means 
that the freeze fault is detected instantaneously. 

The results for several different freeze fault amplitudes on hı and the respective fault detec- 
tion times Age; are given in Tab. 8.1. All detected faults were checked in detail to identify 
settings showing the CMP effect. 

In case of fr = 6.2cm a CMP condition occurred. This is due to the fact that the used 
freeze fault intensity is enclosed in the 6; = 0.05cm interval around the real value of 
his0 = 6.22cm. Thus at k = 651 the freeze fault case is really close to the correct value 
of hı 650+1, meaning that there is a feasible parameter mapping hı, k = 6.2 to hi k41 = 6.2. 
Subsequent points do not provide additional information, as the sensor value is fixed by the 
freeze fault. The only additional information is provided by the pump signal vı. The varying 
pump signal leads to a movement of the center of the outer enclosing zonotope. At kget = 684 
the center of the zonotope is moved to a position that generates a CMP effect that is erro- 
neously interpreted as an inconsistency. 
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Figure 8.12: Verification results for freeze fault of f = 6.0cm at kerr = 650 
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All calculation times for freeze faults given in Table 8.1 are less than the time of the measure- 
ment signal Teate < 2000s = T 


Table 8.1: Different fault amplitudes and resulting detection times for freeze fault 


Fault f; Time step of Quality Calculation Time 
incm Fault kerr | Detection kdet Teale in S 
7.0 650 650 no CMP 1181.0 
6.5 650 650 no CMP 1200.9 
6.2 650 684 CMP occurred 743.4 
6.0 650 650 no CMP 1461.7 


The second regarded malfunction is an offset fault. In this case the sensor value is not fixed, 
but a specific value is added to each measurement: 


Sk Vk € [1, kn = 1] 
ok = (8.41) 
Sk + fo Vk € [Kerr, 1] . 
Again the measurement values of hı and hg are enclosed using ô}, = 67, = 0.05cm. The 


results for an offset fault of f; = 0.7cm on sensor hı are depicted i in Fig. 8.13. Again the 
verification result changes to infeasible right at the moment the fault gets effective i.e., at 
kdet = 650 = kerr- 

Further results for different fault amplitudes are given in Table 8.2. Instantaneous detection is 
possible up to a fault amplitude of fo = 0.15cm. Note that the measurement noise is enclosed 
using ô}, = ô, = 0.05cm which leads to an interval width of 207, = 0.1cm. This is very 
close to ‘the fault amplitude f, = 0.15cm. When using f, = 0. lem - which is exactly the 
interval width - the CMP effect occurs. The necessary calculation time is less than the signal 
duration for all regarded offset fault intensities. 


Table 8.2: Different fault amplitudes and resulting detection times for offset fault 


Fault fo Time step of Quality Calculation Time 
incm Fault kerr | Detection kaet Teale ins 
0.70 650 650 no CMP 1367.9 
0.30 650 650 no CMP 1475.0 
0.20 650 650 no CMP 1361.5 
0.15 650 650 no CMP 1210.0 
0.10 650 650 CMP occurred 1036.8 
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Figure 8.13: Verification results for offset fault of fo = 0.7cm at kerr = 650 
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8.3.3 Multiplicative Faults 


Multiplicative faults can be related to faults in system components i.e. a congested or leaking 
pipe or decreasing pump performance. Such a multiplicative fault fg directly influences the 
system parameter: 


(8.42) 


ENT Vk € [1, kerr — 1] 
"oS Aa qe. Vie [kern T]. 


A maximum absolute deviation of ôp, = op, = 0.05cm is used to enclose the measurement 
values of hı and hg. An exemplary setting for fg). = 0.035 at kerr = 1200 is depicted in 
Fig. 8.14, further results are given in Tab. 8.3. 

It can be seen, that the faulty parameter influences the value of hi. This change in system 
dynamic is recognized by the verification method at kacı = 1200 = kerr. 

All detected inconsistencies were checked in detail. Reliable results are possible up to 
fe, — 0.010. This is a very small value with respect to the nominal parameter variability 
ON = 0.025 which shows the new method is very sensitive. 

The condition Teate < 2000s = T holds for all entries in Tab. 8.3. 


Table 8.3: Different fault amplitudes and resulting detection times for parameter fault 
Fault fg) Time step of Quality Calculation Time 
Fault kerr | Detection kaeı Teale ins 
0.035 1200 1200 no CMP 1036.6 
0.022 1200 1215 no CMP 1123.9 
0.020 1200 1227 no CMP 1166.7 
0.010 1200 1572 no CMP 764.7 
0.005 1200 1638 CMP occurred 746.7 
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Figure 8.14: Verification results for multiplicative fault of fec 1) = 0.035 at kerr = 1200 
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8.4 Application: Diagnosis By Kaucher Based Guaranteed 
Verification (Single-Tank) 


The diagnosis method is now applied to real measurement data instead of simulation data 
as in the previous chapter. Therefore the IRS three-tank setting (introduced in Section 8.1) 
is used, again reduced to the single-tank setup. The respective geometric parameters can be 
taken from Appendix G, Tab. G.1. 

The following scenario is regarded: The water level in tank 2 has an initial height of 
ho = 24.48cm and is rising due to the input flow from pump 2. Pump 2 is running at a 
high load with varying intensity. 

First, the fault free setting is evaluated to show that the method is able to verify the nominal 
setting. Then the additive sensor faults “freeze” and “offset” are applied to the measurement 
data. Finally a scaling fault on the height measurement data is considered. 


The results of different fault intensities as well as the detection and calculation times are 
given in several tables. Each result was evaluated carefully to determine CMP conditions. 
The results were initially published and presented in [Sch18c]. 


8.4.1 Fault Free Setting 


The regarded operation range is defined to be ha € [24, 46] cm and used with (8.5) to obtain 
the nominal range 0* = [0.971, 0.979]. 

The result is calculated by using an interval width of ô, = 0.4cm to enclose the measure- 
ment data of ha. The pump measurement data is assumed to be noiseless and thus used as 
point real value. 

The fault free behavior is depicted in Fig. 8.15. It is verified for the entire measurement time. 
The necessary calculation time is Teale = 34.4s < 120s = T. 
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Figure 8.15: Fault free measurement data of the single-tank diagnosis scenario 
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8.4.2 Additive Faults 


The first considered additive fault is given by a sensor freeze. The used mathematical model to 
distort the fault free measurement data of ha is given by (8.40). The measurement is enclosed 
using an absolute deviation of 07, = 0.4cm. The resulting system run for fr = 37.9cm on 
the measurement of ha at kerr = 60 is depicted in Fig. 8.16. 

It can be seen that it is not possible to verify the measurement data as soon as the freeze 
fault is active. The failure is detected at kget = 60 = kerr, i.e. at the very first time the 
measurement is distorted. 

An evaluation of the performance of the method for several different freeze fault intensities 
Fr is listed in Tab. 8.4. It can be seen that it is possible to detect faults in a large range 
from f; = 42.0cm to ff = 37.9cm. The lower value is very close to the correct value 
he hep. = 31.54cm. 

All faults are detected directly at their first appearance, i.e. at kget = 60 = kerr. All results 
were checked carefully to ensure that there is no CMP effect present in the results. 

All calculation times are less than the measurement time, i.e. Teale < 120s = T. 


Table 8.4: Different freeze fault amplitudes for sk., = h2,60 = 37.54cm 


Fault f; Time step of Quality | Calculation Time 
incm Fault kerr | Detection kj, Teale in S 
42.0 60 60 no CMP 99.6 
39.0 60 60 no CMP 49.6 
38.5 60 60 no CMP 50.2 
38.0 60 60 no CMP 80.8 
37.9 60 60 no CMP 31.0 
37.7 60 not detected | no CMP 24.9 
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Figure 8.16: Measurement data with freeze fault fr = 37.9cm 
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Second, an offset fault setting is applied. Therefore a constant offset f, is added to the faultless 
measurement data of ha, according to (8.41). Again, an absolute deviation of 07, = 0.4cm 
is used to enclose the measurement. The performance of the method is shown exemplary in 
Fig. 8.17 and Fig. 8.18. 

The large offset of f, = 5cm in Fig. 8.17 is rather obvious and could also be detected by an 
expert. On the other hand, the very small offset of f, = 0.35cm in Fig. 8.18 is very hard to 
distinguish from the fault free measurement depicted in green. 

Nevertheless the zonotopic method is able to detect it at the moment of its first appearance. 
This is a very powerful property as the detected offset of f, = 0.35cm is smaller than the 
used interval radius of 57), = 0.4cm. 

This performance is due to the dynamic between time instant kerr — 1 and time instant kerr 
that is created by the appearance of the offset fault. This dynamic is detected instantaneously 
as it is outside of the nominal parameter range. 


Several results for different offset intensities are given in Tab. 8.5. The table shows that in- 
stantaneous detection, ie. kaet = kerr is possible within the range of f, = 0.35cm and 
fo = 5cm. No result shows the CMP effect which means that they are of good quality and 
provide reliable results using a zonotopic approximation of the united solution set. 

The calculation time of all results is also given in Tab. 8.5. It can be seen that the measure- 
ment data can be processed in less than the genuine signal time, i.e. Teate < 120s = T holds 
for all fault intensities. 


Table 8.5: Different offset fault amplitudes for sp = h2 60 = 37.54cm 


err 


Fault f, Time step of Quality | Calculation Time 
incm Fault kerr | Detection kae; Teate in S 
5.00 60 60 no CMP 97.2 
2.00 60 60 no CMP 92.0 
1.00 60 60 no CMP 59.3 
0.50 60 60 no CMP 60.1 
0.35 60 60 no CMP 57.0 
0.20 60 not detected | no CMP 25.2 
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Figure 8.17: Measurement data with offset fault fo = 5cm 
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Figure 8.18: Measurement data with offset fault fo = 0.35cm 
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8.4.3 Scaling Faults 


In the diagnosis scenario discussed in this section so far, real measurement data from a single- 
tank process is used. To realize multiplicative faults with the same measurement data, a 
scaling fault in the corresponding sensor of ha is assumed. This leads to the following scaling 
fault model 


Ssk = TEE ee M (8.43) 


Sk ' fs Vk = [Kerr, T] 


which replaces the former multiplicative fault (8.42). The absolute deviation to enclose the 
measurement of ha remains ô}, = 0.4cm. The results are depicted in Fig. 8.19 and Fig. 8.20 
for f, = 0.95 and f, = 1.01 respectively. 

It can be seen that even factors very close to one (e.g. fs = 1.01, meaning a deviation of 1%) 
can be detected. 


Results for an extensive range of factors are given in Tab. 8.6. 

It is not possible to detect the fault intensity of fs = 0.97 as this parameter is very close to 
the nominal parameter 0* = [0.971, 0.979] representing the desired system dynamics. This 
means there is a deviation of 0.1% between f, = 0.97 and 0* = 0.971 which is one order of 
magnitude less than for f, — 1.01. 


All successfully detected faults lead to kaget = 60 = kerr, i.e. they are detected right at their 
appearance. There was no CMP condition present in the regarded settings. 

Again, it is possible to calculate the results for all fault amplitudes in less than the genuine 
signal time, i.e. Teale < 120s = T. 


Table 8.6: Different scaling fault amplitudes for s7.,,.,. = h2,60 = 37.54cm 


Fault fs Time step of Quality | Calculation Time 
Fault kerr | Detection kae; Teate in S 
1.10 60 60 no CMP 77.7 
1.05 60 60 no CMP 80.4 
1.03 60 60 no CMP 77.3 
1.01 60 60 no CMP 64.6 
0.97 60 not detected | no CMP 31.9 
0.95 60 60 no CMP 108.8 
0.90 60 60 no CMP 96.6 
0.75 60 60 no CMP 75.2 
0.50 60 60 no CMP 117.8 
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Figure 8.19: Measurement data with scaling fault f; — 0.95 
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Figure 8.20: Measurement data with scaling fault f; — 1.01 


8.5 Conclusion 143 


8.5 Conclusion 


The methods and theories developed throughout this thesis were applied and demonstrated 
in this chapter. 


First, simulation data of a single-tank process was used to show the performance of the ver- 
ification method based on Kaucher arithmetic and zonotopic inner enclosures of the united 
solution set. 

This approach was extended to the hybrid setting given by a two-tank system. It was shown 
that the introduced method is able to verify the correct system in case of known switching 
times and known active subsystems. 

The application to a four-tank process showed that it is not possible to verify the system in 
various faulty settings even for very small fault amplitudes. This is a relevant indicator that 
a fault is present in the system and can thus be used for fault detection. The performance of 
the developed method was shown for three different and common fault types, namely freeze 
fault, offset fault and multiplicative fault. 

The same diagnosis algorithm was finally applied to real measurement data provided again 
by the single-tank process. It could be shown that the algorithm obtains valuable results by 
detecting even very small faults from real world data. 


The calculation time of all introduced examples was less than the genuine time of experi- 
ment on a standard laptop. Therefore the application might in general be suitable for online 
application in a diagnosis setting. 


9 Conclusion 


Modern engineering is able to develop and build complex and powerful systems to an un- 
precedented extend. The functionality of these systems is rapidly increasing and masters 
tasks that used to be subject to highly trained humans. The challenge how to build such sys- 
tems is nearly completed. Still remaining is the question how to ensure correct functionality 
of such powerful safety critical systems. Current safety analysis relies on sophisticated meth- 
ods from the field of testing. Even though these methods are very mature, they are essentially 
falsification approaches meaning that there are type II errors by definition. However, in the 
case of safety critical systems, it is necessary to ensure the absence of type II errors. 

This thesis provides the foundations for a new specification and verification approach able 
to provide the necessary type II error free results. 


Therefore a new notion of set based consistency for dynamic systems with a given speci- 
fication is presented. Kaucher interval arithmetic is used to enclose the measurement data 
in a bounded error sense. Thus, the specified behavior of a dynamic system can be verified 
by measurement data even in the presence of noise and sensor uncertainty. Consistency is 
defined using the Kaucher arithmetic united solution set which leads to mathematically guar- 
anteed results. The verdicts calculated by the new Kaucher based method can not show type II 
errors (hidden faults) by definition and are thus suitable to provide a reliable verification of 
safety critical systems. 


It was proven mathematically that this holds for a wide class of systems, including time in- 
variant, interval type and hybrid systems, which can be used to describe even nonlinearities. 
The notion of consistency was extended to include the discrete event part of a hybrid system 
and requirements on the connection of the two system classes were derived. Several exten- 
sions were introduced, leading to a new iterative identification and segmentation algorithm 
for hybrid systems which is able to handle even unknown switching times. In case the calcu- 
lations can be done fast enough, the developed approach can also be used for the diagnosis 
of dynamic systems. Requirements on sampling time and hardware performance have to be 
determined for each specific setting individually. 


The presented methods were successfully applied to several example systems, consisting of 
a variation of different tank settings. The results were shown, interpreted and discussed. 


The results provide the base to answer the research question that governed this thesis. The 
new theories, methods and algorithms developed in this thesis form the foundation for reli- 
able safety analysis of highly automated safety critical systems. The results of this thesis can 
be used to solve the arising problems of current powerful and interconnected systems that 
are increasingly interleaving our daily live. 


A Analysis Perspectives 


A popular definition to distinguish different analysis perspectives was coined by [Boe84]: val- 
idation means “building the right product” whereas verification means “building the product 
right”. 

This is used to set up a high level differentiation between validation and verification / falsifi- 
cation as given in Fig. A.1. Validation is always concerned with the desire of the customer and 
evaluates the question whether the developed functionality fulfills this desire. The field of val- 
idation is very important and large research effort including psychology, behavioral science 
and linguistics has been put on it in the last decades [Mac95][Que98][Fau03][Fol08][Bar13]. 
Nevertheless validation is not part of this thesis. 

The objective of security focuses on the detection of intentional misuse by (un)authorized 
subjects, e.g. due to hacking attacks or user errors. The whole field of security, including 
conscious misuse, hacking or manipulation is also not in the scope of this thesis. 


«— — Development Process ————~>]<— Runtime —> 


Customer Validation 


Specification : 
H id Falsification 
uman 


: Implementation 


Diagnosis 


Figure A.1: Evaluation terminology 


The realization of the specification is done by human engineers, therefore it is likely that 
there are mistakes during the process of implementation. A wide spread approach to find 
these mistakes is given by the concept of falsification, which tries to determine a so called 
counter example that shows unspecified or wrong behavior. If it is not possible to determine 
a counter example it is assumed that there are no counter examples at all and thus the system 
is considered to be correct. However, due to limited runtime of the falsification process it is 
possible that there are undetected (hidden) faults in the implementation. Therefore type II 
errors are possible which is an disadvantage in case of safety critical systems. 

If the specification is very formal, the implemented system can be analyzed in a formal way. 
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Verification methods aim on proving the correctness of the implemented system in (all) op- 
erating conditions. The goal is to prove that the system always shows nominal behavior. 
Verification and falsification methods are in general conducted during the development pro- 
cess, while the system operates in some kind of artificial environment. The evaluation can be 
done offline and might thus need more calculation time or can be run several times during the 
development process. Mistakes occurring during system operation are tackled by methods 
of the diagnosis and monitoring field. They need to run online in parallel to the real system 
operation and are thus required to be very fast. 

In case of model based diagnosis, a model of the nominal system is generated that is used to 
calculate the nominal system output in parallel to the real verification object (VO). Therefore 
the inputs of the VO are measured and also applied to the nominal model. The resulting out- 
puts are compared with the measured outputs of the VO which leads to a so called residual 
vector (see among others [Ise93][Ven03a] [Ven03c][Ble10]). In case of an undisturbed system, 
the residual vector is zero if there is no fault present in the VO. A fault is detected if there 
is a non-zero residual vector. If it is necessary to gain further knowledge of the fault, more 
sophisticated methods can be applied to localize the exact point of fault occurrence within 
the VO ([Ven03a][Ven03c][Che14]). A drawback of the diagnosis approach is that - due to 
measurement noise and model imprecision - the residual vector is not always exactly zero 
even in the fault free case [Ise06, p. 198]. 


B  Derivation of the Interval Distribution 


This appendix provides the derivation of a probability distribution on the parameter p con- 
necting two intervals u and y. 

Two examples are presented in Chapter 3. Example 3.4 shows a setting with a proper result 
of p and Example 3.5 demonstrates a setup leading to an improper solution. The interval 
ranges of u and y are sampled with Au = Ay = 0.0001 and used to calculate the resulting 
parameter p, for all possible combinations. 

It is also possible to theoretically derive the shown results. Therefore, two random variables 
u and y are defined with uniform distribution between the infimum and the supremum ofthe 
interval values u and y. The probability density functions of the two random variables are 
given by: 


piel u-u tees (B.1) 


-y == (B.2) 


A general probability density function according to [Bro08, p. 816] has to fulfill the assump- 
tions 


f (x) 2 0, Va (B.3) 
L. f (x) da — 1. (B.4) 


Assumption (B.3) is valid for (B.1) and (B.2) by definition. Assumption (B.4) can be shown as 
follows: 


zd fa (u) du4- fuo) dus [ fa (u)du 


Sa eS ON 
0 
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The proportional parameter p with 
up=y (B.5) 


thus can be interpreted as random variable 


y 
p — g(u,y) — —. (B.6) 
u 
There are different possible realizations of u and y depending on the specific values of p. 
Therefore the probability density function of p is according to [Jon02, p. 118] given as 


fi = f ee (5.7) 


With the constant densities of the uniform distributions f, (u) and f, (y), and assuming only 
non-negative input values u > 0, (B.7) can be relaxed to 


oo 
RO= [u fale) fuas Bs 
0 — ——— 
constant for constant for 
u<usu YJu<p<Ylu 
else O else O 


The antiderivative is zero for all p ¢ [v/u, ¥/u] and !/p ¢ [v/y, T/y]. Else the densities consist 
of constant values, leading to the antiderivative 


1 oo 
tp (p) = sese? (B.9) 
0 
with c, — — and c, = = The evaluation of f, (p) depends on the infimum and supre- 


mum of wand y and can be generalized as 


fo (p) = 366r max (o (min(e, maxi). e c min (s/p, v9) | 
(B.10) 


It is now possible to draw the derived density function for specific values of u and y. Exem- 
plary plots for a proper and an improper setting are given in Example B.1. 
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Example B.1: 

The plot for the values of u = |2, 3] and y = [4, 9] according to Example 3.4 is depicted 
in Fig. B.1. The density function has the same shape as the sampling based result given in 
Fig. 3.5. The plateau in the figure shows that the solution is proper. 
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Figure B.1: Probability density function fj (p) for the proper case 
The plot of an improper setting according to Example 3.5 is given in Fig. B.2. The input and 
output intervals are u = |2, 3] and y = |4, 5] leading to an improper solution and an eroded 
plateau. 
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Figure B.2: Probability density function f, (p) for the improper case 
These results are in accordance with the theory introduced in Chapter 3. 


C Full Rank Criteria 


It is in general N P-hard to determine whether a given interval matrix A has full rank, re- 
spectively to check the matrix for singularity [Sha14]. However, there are some criteria to 
determine the property of full rank [Sha14]. Four sufficient conditions and one necessary and 
sufficient condition are given in the following. 

It is necessary to introduce the absolute value of an interval 


|z| = max (|| , æl) (C.1) 
and the magnitude 
"- min (|z|, |z|) ,if0 da (C.2) 
0 , else. 


The first sufficient condition for quadratic problems is based on diagonal dominance. The 
interval matrix A € IR'*") is nonsingular, if it is diagonal dominant. This means the 
inequality 
n 
a+ >) | (C.3) 

j=l 

i#i 
holds for i € {1,2,...,n}. 
There are two approaches to extend this condition to overdetermined equation systems, i.e. 
A € IR(”*”) with m > n. The first approach searches for diagonal dominant subsquares 
within the overdetermined interval matrix. If there is such a diagonal dominant subsquare, 
the whole interval matrix has full rank. However there might be no diagonal dominant sub- 
square even though the matrix has full rank. This can be due to permutation of rows of the 
matrix. Even though permuted lines do not change the rank of a matrix, it does change the 
appearance of diagonal dominant subsquares. However, the property that permuting rows 
does not change the rank of the matrix can also be used to solve the problem. The lines 


can be permuted algorithmically such that diagonal dominant subsquares are created. This 
condition is still sufficient for overdetermined systems. 
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A second sufficient condition for full rank of an interval matrix A € IR(”*”) is based on 
the spectral radius. Thereby the spectral radius p (A) is defined to be the largest absolute 
singular value of the matrix A [Lax02, p. 195]. If the spectral radius fulfills 


p(|(4e)"] Aa) «1 (C4) 


and the center matrix A, has full rank, also the interval matrix A has full rank. Thereby 
A= (ATA). AT denotes the pseudo inverse of the matrix A € R(?*") with m > n. 


A third sufficient condition is based on the singular values of the matrix A € IR”), If the 
condition 


Omax (AA) X. O min (A.) (C.5) 


is fulfilled, the interval matrix A has full rank. Thereby omax (A) and omin (A) denote the 
greatest, respectively smallest, singular value of the matrix A. The singular values are defined 
as the nonnegative solutions to the system 


(a 4) (eG). es 


Condition four uses an absolute subordinate matrix norm || - || of A € IR*"). Assuming 
full rank of the center matrix A., the sufficient condition is given by 


[|Aal] < [AXI (C.7) 


If (C.7) holds, A has full rank. The proofs of all four sufficient conditions are given in 
[Sha14]. 


According to [Roh12] there is a fifth, necessary and sufficient condition: An interval matrix 


A € IR*” with m > n has full rank iff 
|A-X| € Aa |X| (C.8) 


with X € R” can only be solved by the zero solution X = [0,0,...,0]". 

Sufficiency is based on the idea that there is a non trivial solution X ¥ [0,0,... u as 
soon as the matrix A does not have full rank. Necessity follows from the existence of the 
non trivial solution. If this is the case, A cannot have full rank or the non trivial solution 
does not solve (C.8). An extensive proof of this condition is given in [Sha14]. However, the 
approach directly aims on an N P-Hard problem which means that it can only be checked 
approximately. 


D Existence and Uniqueness of the Algebraic 
Solution Set 


The sufficient conditions on the existence of an algebraic solution given in [Sha96], [Mar99] 
and [Lak99] are sketched in this appendix. To follow those ideas, two more interval arithmetic 
notations are necessary. 

Using the dual (-) operator given in (3.37), the proper projection pro (a) is defined as 


fri 
stu x ‚if is proper (D.1) 
dual (x) , else. 


The second property is ı-nonsingularity. A quadratic point real matrix Q € R(**") is called 
ı-nonsingular if 


Qzx =062=0€ IR” (D.2) 


holds. Otherwise Q is called »-singular. 


According to [Sha96] there is an algebraic solution 57, to the interval linear equation 


Ar = B with A € IR™*” for any B € IIR", if A is sufficiently narrow and pro (A) 


contains an ?-nonsingular point matrix. 
Thereby “sufficiently narrow” means that | AA | is sufficiently small. 
The proof of existence given in [Mar99] is based on the iterative approach to determine the 


algebraic solution set given in [Kup95]. For this proof, the notation of the diagonal matrix 
D(A) is introduced for an interval matrix A € IR®*™ 


a3 .ifi-j 


= (UNO (D.3) 
1<i<nl<j<n 0 fd X 3. 


D(A) — (a?) 
The inverse of the diagonal matrix is given by 


E = Fij) = l/dual(a (2) 5 ifi= j 
E (4) - (d P f ies j. en 


Stayt Mp mm 
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Using the dual (-) operator from (3.37). The iterative solution algorithm given in [Kup95] 
converges to the algebraic solution 5 7, if 


ID (ls (D.5) 
IA + opp (D(A)) || <1 (D.6) 
holds, with opp (-) according to (3.34). The used matrix norm || - || is the maximum of the 


linewise sum of the absolute interval values (C.1): 


All — ex (>: ja") | = 
> k=1 


The interested reader is referred to [Mar99] for further considerations. 


A generalized approach for overdetermined systems A € IR(™*") was introduced in [Lak99]. 


The regressor matrix A is split in three parts with A = Ap + Aı + Ag and 


G4) > 0 
a= (al) ana? D8) 
<i<m,1<j<n 0 E else 
T G3) ifai) eg 
DC NNNM EM MEE (D.9) 
<i<m,1<j<n 0 3 else 
- G3) Gf glia) qi) 
Bela. |e dida (D.10) 
<i<m,1<j<n 0 E else. 


The problem can be reformulated as an extended system that considers the upper and lower 
bounds of the interval values explicitly, as given in [Lak99]. This problem can then be trans- 
ferred to a set of inequality conditions. It is possible to show that there is not more than one 
solution for any b € B ifthe derived set of inequality conditions has zero as unique solution. 
The extensive proof is given in [Lak99]. 


E System Behavior Specification 


The verification methods developed in this thesis are based on the assumption of a system 
specification available in ARX form. In a practical setting, it is necessary to determine these 
nominal parameters. Control engineering specifications are in general based on tolerance 
bands, steady-state errors, rise and settling times or acceptable overshoots. Such specifica- 
tions inherently show interval properties - even though in general no interval arithmetic is 
used. 

This appendix provides two approaches to determine the ARX parameters of such intuitive 
graphic specifications. In the time domain, a method is introduced to determine the parame- 
ters from a desired step response. This step response can be set up using the drag-and-drop 
function provided by a toolbox. A second method determines the parameters from the fre- 
quency domain. Therefore only the tolerance band widths and pass/cut-off frequencies of a 
filter function need to be specified. In case the Kaucher based method is applied in a diagnosis 
setting, it is beneficial to use a nominal physical model of the regarded process. This physical 
model can also be used to determine the desired ARX parameters. 


E.1 Time Domain Specification 


This specification approach allows an intuitive specification of the desired behavior based on 
time domain input-output behavior. The first step is to specify an input signal and the desired 
resulting output signal. Also the class of the desired system behavior has to be given. The time 
domain input-output behavior is then used to determine the parameters of a transfer function 
in the complex s-plane. In this appendix a step input is used to determine the properties of 
a proportional gain first order time delay system (PT1). The method given in [Föl13, p. 77] 
can be used to determine the system parameters based on a given step response. A time 
continuous step response denotes the output values y(t) generated by a step input 


0, k«0 
t)= ot) = E.1 
u(t) = e(t) n ie (E.) 
The complex frequency domain transfer function of a basic PT1 system is given by 
k 
G(s) — - E.2 
(S) 1+Ts e» 


with gain k, and time delay T. If both parameters are determined, G(s) can be transformed 
to its discrete time representation. The desired nominal parameters ©* are then given by the 
parameters of the discrete time transfer function. 
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For given input-output data [u(t), y(t)] with t € [0, T] itis possible to determine the transfer 
function parameters. The gain k, is given by the stationary value k, = Yoo which is defined 
to be the last point of the measurement y(T), assuming T is large enough to allow y(t) to 
settle. If the output signal is sampled with sampling time At, the information of the resulting 
points can be used to calculate the time delay T. Therefore each available sampling point 


Uk = y(kAt) is used to calculate an auxiliary value 


me =1- 2 (E.3) 
Yoo 
which is then used to determine the time delay 
~ kAt 
Tk =— E.4 
= m) = 
The time delay of the transfer function T is given by the arithmetic mean of all nj, = E 
values of Tus ie. 
> WA 


Afterwards the transformation to discrete time is done. The resulting parameters O* are used 
to determine the ARX step response. 

The introduced functionality is implemented in a Toolbox. The application of this Toolbox is 
demonstrated in Example E.1. 
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Example E.1: 


Based on an initial arbitrary input-output signal, the toolbox provides the possibility to 
move the sampling points via drag-and-drop. Fig. E.1 shows an exemplary output signal 
that was created as the step response of a PT1 system with gain kp = 10 andT = 2. 
The resulting continuous trajectory was sampled with At = 1s. The blue points depict the 
sampling points that can be moved using drag-and-drop. The user can move the sampling 
points such that the resulting trajectory shows the desired behavior. In this case, the toolbox 
calculates the system parameters according to (E.3)-(E.5). The resulting system is able to 
generate the desired values for the given step input. The depicted setting leads to the time 
discrete ARX parameters a = 0.64227, c = 3.6077. The respective time discrete step response 
is given by the red trajectory in Fig. E.1. 


10 - 


Output signal 
Or 
T 


— Specification 


— ARX 
0 tc | | l | | l | n L 1] 


0 1 2 3 4 5 6 7 8 9 10 
Time step k with AT = 1s 


Figure E.1: Time domain specification toolbox 
It can be seen in Fig. E.1 that the time discrete trajectory (red) is close to the specified points 
(blue). This demonstrates that it is possible to determine the ARX parameters of a system by 
using a graphical user interface with drag-and-drop to set up the desired step response. 


E.2 Frequency Domain Specification 


This specification approach allows an intuitive specification of the desired behavior based 
on designing the amplitude response of the system. First order systems can be interpreted 
as low pass filters. Each Filter has a specific frequency domain characteristic consisting of 
the location and the width of the passband and the stopband. Based on this information, the 
method of [Lüc80, p. 147ff] is used to determine the filter coefficients. The resulting filter can 
be transformed to discrete time which leads to the desired ARX coefficients. 
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The method is applicable for low-pass, high-pass, bandpass and band-rejection filters. In 
this appendix the design of a low-pass is presented. Therefore the cut off frequency Qp, 
the stop band frequency ©), and the respective passband width A, and stopband width A, 
(see Fig. E.2) need to be defined by the user. These frequencies are defined with respect to the 
periodic interval of the frequency response 0 < €) € 7. The frequencies are now transformed 
into the 0 < f < co domain by 


fp = tan (®»/2) (E.6) 

fs = tan (8/2) . (E.7) 
These values lead to the normalized low-pass representation 

Tonart =1 (E.8) 

Fs;norm = fsf/ fp. (E.9) 


The normalized low-pass representation can be achieved for all four kinds of filters, by using 
different transformations. The following design routine is thus applicable in every setting. 
To ensure that the given passband and stopband limits are met, the auxiliary variables 


E 2A4— AZ 

Aq = V Ti for 0 < f <1 (Passband) (E.10) 
1— Aq 

xo cL AD 

As = EN EE for Tenor < f (Stopband) (E.11) 


are calculated. The transfer function of an exponential filter is then given by 
Ch = Goma r r (E.12) 


with poles f,.,; and normalization constant Go. The denominator order q € N can be calcu- 
lated with 


> logio (-/A,) 


> : (E.13) 
log10 (fs norm) 
The real part and the imaginary part of a complex pole fx ; are given by 

E 2i —1 
X (fai) = —€ sin ( : ; z) (E.14) 

z 2i—1 
S(fooi) = € V4 cos ( f z) ] (E.15) 

q 


The factor e and the respective normalization can be used to adjust the level of the frequency 
response. The used value of e can be chosen from the interval 


A, x 


Ganon? ^’ en 


e= [e g = 
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Thereby choosing c = € means that the frequency response touches the constrained region at 
the end of the passband fp norm, whereas e = € means that it touches the constrained region 
at the beginning of the stopband fs norm- 

The normed lowpass is now transformed back to its genuine frequency form and afterwards 
into the time discrete representation in order to extract the desired nominal parameters O* 
of the transfer function. An application of this method is given in Example E.2. 


Example E.2: 

The frequency domain specification of a low-pass filter is given by A, = 0.1, A, = 0.2, 
N, = 0.4» and), = 0.77. This means that the desired frequency response is located within 
the green area in Fig. E.2. The introduced filter design procedure of [Lüc80, p. 147ff] is applied 
to the setting. 

Based on the specified values, the choice e = € leads to the ARX coefficients 


[a1, a2, a3] = [0.1425, —0.3387, 0.0130] (E.17) 
[c1, €2, €3, Ca] = [0.1479, 0.4437, 0.4437, 0.1479] . (E.18) 


The respective frequency response is depicted as solid blue line in Fig. E.2. 
Using the same values but choosing € = e leads to 


[a1, a2, a3] = [-0.2643, —0.3518, —0.0244] (E.19) 
[c1, c2, c3, c4] = (0.2051, 0.6152, 0.6152, 0.2051] , (E.20) 


displayed as dashed line in Fig. E.2. 
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Figure E.2: Frequency domain specification toolbox 


F Excitation Signal Design 


Hybrid system verification poses specific requirements on the excitation signal. It is assumed 
that these requirements are fulfilled throughout this thesis. However, determining a suitable 
excitation signal is in general not trivial. The excitation signal of a measurement is often 
chosen depending on the intended purpose of the experiment. Arbitrary noise signals (white 
Gaussian noise) can be used to ensure persistent excitation of all frequencies. Arbitrary mean- 
ingful signals (impulse or step signal) are used to perform control theoretic modeling such as 
impulse response or step response. Also there are specifically designed input signals fitted to 
the implemented logic in the current verification object. 

In the context of hybrid systems as regarded in this thesis, there are two properties that need 
to be fulfilled. Each subsystem with its respective individual dynamic needs to be persistently 
excited. Furthermore, all states of the superimposed state machine need to be activated once. 
Therefore the respective switching thresholds have to be met to enable the switch event. The 
situation that a switch is triggered during the excitation and identification phase of each sub- 
system has to be avoided. A first possible solution idea was developed in the master thesis 
[Rie17]. The basic outline is sketched in this appendix. The method uses three steps: 


1. Path calculations to ensure state coverage of the superimposed state machine 
2. Design of a persistent excitation signal without leaving the subsystem 


3. Efficient transfer of the subsystem to its switch threshold 


F.1 Path Calculation 


The superimposed state machine is transformed to its graph representation Gz. A coverage 
algorithm is used to determine paths that include all states and all transitions of the graph. 
Additionally, the length of the path needs to be minimal to enable short measurement times. 
If such an optimal excitation signal is used, missing states or transitions can be used to prove 
inconsistency. 

The problem of state and transition coverage can be reduced to transition coverage only. 
This is due to the structure of the specification, where each state needs to be connected to a 
transition. 


XL F Excitation Signal Design 


To use a modified depth first state coverage algorithm, the graph Gz is transformed such 
that all transitions are represented by states in the transformed graph G’ z and vice versa (see 
Fig. F.1). 


XA 
© 


Figure F.1: Genuine graph Gz (left) and transformed graph G’ z (right) 


A modified recursive depth first algorithm is started in a specified initial state and checks the 
number of possible successors of each successor of the current state. Also the distance to the 
successors are taken into account to enable short paths. The result of the algorithm is a path 
that covers all states of Gy and thus all transitions of Gz. This result is used to identify each 
subsystem in the path. 


F.2 Persistent Excitation Based on Fisher Information 
Matrix 


Persistent excitation of each subsystem is ensured by a specific input signal. This input signal 
is calculated based on the Fisher information matrix [Eba14][Man10], which is only applicable 
for stable systems. Using the Fisher information matrix M, the parameter covariance of an 
estimator is limited by the Cramer Rao Bound [Goo77] to 


cov(©) > M (F.1) 


F.3 Transfer to the Switch Threshold XLI 


Thereby M is defined using the expectation E {-} as 


8log (p (y | a [2s (p(y | a) 


06 86 c 


M(0,U)=E | 
The probability log p(y | O, U) resembles the situation that y is observed if the true parame- 
ters are given by © while using the input U = (ux) ae To achieve a parameter covariance as 
close as possible to the Cramer Rao Bound, the Fisher matrix has to be maximal with respect 
to the input signal used and the parameters. This can be achieved by using D-optimality for 
the given nominal parameters as defined in [Man10]: 


U* = —min(log det(M(0*,U))). (F.3) 


To solve the optimization Problem F.3, an initial feasible input signal U;,;; is chosen. This 
signal is then optimized iteratively for each time step uj; ;; ;, until the optimization converges. 
Each input value u% is thereby bounded to the range of feasible input values given by the 
user. 


F.3 Transfer to the Switch Threshold 


After the identification of the subsystem, it is necessary to activate the successive switch. 
This is done by transferring the relevant system value within its activation limits 1. 

The specific event and thus the activated transition are already determined in the result 
of the path calculation. This is done using the well known Hamilton formalism. There- 
fore the objective function is set up in terms of the difference between the desired value 


yo = 5 & + i" and the current value yx, i.e. AYk = Yk — Y@- The resulting objective 
function is given by 


T-1 


1 1 
J =5AyrSAyr +5 2. Ayr Q Ayr (F.4) 


with the penalty matrices S = — and Q = 1. It is now possible to set up and solve the 
Hamilton equations [Sag68]. The first step is to transfer the ARX system description to a 
vector matrix notation: 


Yk+1 ay ag ... ang Uk Ci ... Cn. Uk 
Uk 1 0 ... 0 Uk—1 Ü c 0 Uk—1 
; —-|. . . . tl. o. . . (E5) 
Uk—na-4-2 0 ... 1 0 Yk—natl 0 ... 0 Uk—n.t+l 
| 
A Y. C Uy 


Then the Hamilton equation is set up: 


1 = " 
H (yr, ur, k) = 5A SAy + Akı (AYs + CU;). (F.6) 


XLII 


Excitation Signal Design 


The derivative of H is given by 


OH 
zent. 
Ou 
with 
de = Ayk + AT Agua 
leading to 
H 
Au= -a a 


ult) = ul + Aug. 


(E.7) 


(F.8) 


(F.9) 


(F.10) 


The parameter o is used to scale the result in case the calculated solution violates the feasible 


input range. 


This procedure is applied to all subsystems within the calculated path to construct the overall 
excitation signal. This signal is then applied to the VO and the resulting output values are 
measured. The resulting input and output measurement data can then be used in any of the 


methods introduced in this thesis. 


G Tables of Geometric Parameters 


The parameters of the three-tank lab setting at the Institute of Control Systems (IRS) are 
given in Tab. G.1. 


Table G.1: System properties of the IRS three-tank lab setting 


Value Unit Property 
haau 30.0 cm Height of upper connection valve 
loi 0.0 cm Height of lower connection valve 
ay 0.5 cm? Cross section nominal outflow tank 1 
ag 0.5 cm? Cross section nominal outflow tank 2 
a3 0.5 cm? Cross section nominal outflow tank 3 
A13u 0.5 cm? Cross section upper connection valve vigu 
a13l 0.5 cm? Cross section lower connection valve v13, 
a32u 0.5 cm? Cross section upper connection valve 732, 
4321 0.5 cm? Cross section lower connection valve v35; 
Üleak 0.8 cm? Cross section leakage outflow (only tank 2) 
A1 154.0 cm? Cross section tank 1 
A» 154.0 cm? Cross section tank 2 
A3 154.0 cm? Cross section tank 3 
g 981.0 cm/s? Gravitational force 
yı 16.7 mincm?/(ls) | Constant tank 1 
Ya 16.7 mincm?/(ls) | Constant tank 2 


XLIV G Tables of Geometric Parameters 


The parameters of the four-tank simulation setting are given in Tab. G.2. 


Table G.2: System properties of the simulated four-tank lab setting 


Value Unit | Property 
Vini 0.7 Valve 1 flow to tank 1 
ay 0.071 cm? Cross section nominal outflow tank 1 
ag 0.071 cm? Cross section nominal outflow tank 2 
a3 0.071 | cm? Cross section nominal outflow tank 3 
a4 0.071 cm? Cross section nominal outflow tank 4 
Ay 28.0 cm? Cross section tank 1 
As 28.0 cm? Cross section tank 2 
As 28.0 cm? | Cross section tank 3 
A4 28.0 cm? Cross section tank 4 
g 981.0 cm/s? | Gravitational force 
v1 3.33 | cm/s | Geometric constant 
y2 3.33 | cm/s | Geometric constant 
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